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Using our new numerical-relativity code SACRA, long-term simulations for inspiral and merger of 
black hole (BH)-neutron star (NS) binaries are performed, focusing particularly on gravitational 
waveforms. As the initial conditions, BH-NS binaries in a quasiequilibrium state are prepared in a 
modified version of the moving-puncture approach. The BH is modeled by a nonspinning moving 
puncture and for the NS, a polytropic equation of state with F = 2 and the irrotational velocity field 
are employed. The mass ratio of the BH to the NS, Q = Mbh/Mns, is chosen in the range between 
1.5 and 5. The compactness of the NS, defined by C = GAIj^s/c'^R'ns, is chosen to be between 0.145 
and 0.178. For a large value of Q for which the NS is not tidally disrupted and is simply swallowed 
by the BH, gravitational waves are characterized by inspiral, merger, and ringdown waveforms. In 
this case, the waveforms are qualitatively the same as that from BH-BH binaries. For a sufficiently 
small value of Q < 2, the NS may be tidally disrupted before it is swallowed by the BH. In this case, 
the amplitude of the merger and ringdown waveforms is very low, and thus, gravitational waves are 
characterized by the inspiral waveform and subsequent quick damping. The difference in the merger 
and ringdown waveforms is clearly reflected in the spectrum shape and in the "cut-off" frequency 
above which the spectrum amplitude steeply decreases. When an NS is not tidally disrupted (e.g., 
for Q = 5), kick velocity, induced by asymmetric gravitational wave emission, agrees approximately 
with that derived for the merger of BH-BH binaries, whereas for the case that the tidal disruption 
occurs, the kick velocity is significantly suppressed. 
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I. INTRODUCTION 

The coalescence of black hole (BH)-neutron star 
(NS) binaries is one of the most promising sources for 
kilometer-size laserinterferometric gravitational-wave de- 
tectors such as LIGO [J and VIRGO [5]. A statistical 
study based on the stellar evolution synthesis suggests 
that the detection rate of gravitational waves from BH- 
NS binaries will be about 1/20-1/3 of the rate expected 
for the merger of the NS-NS binaries in the universe 
Then, the detection rate of such systems will be ^ 0.5-50 
events per year for advanced detectors such as advanced- 
LIGO 0, and hence, the detection is expected to be 
achieved in the near future. For clarifying the nature 
of the sources of gravitational waves and for extracting 
their physical information, theoretical templates of grav- 
itational waves are necessary for the data analysis. For 
theoretically computing gravitational waveforms for the 
late inspiral and merger phases of BH-NS binaries, nu- 
merical relativity is the unique approach. 

The final fate of coalescing BH-NS binaries is divided 
into two cases depending primarily on the BH mass. 
When the BH mass is small enough, the companion NS 
will be tidally disrupted before it is swallowed by the BH. 
By contrast, when the BH mass is large enough, the NS 
will be swallowed by the BH without tidal disruption. 
The latest study for the BH-NS binary in a quasiequi- 
librium indicates that the tidal disruption of the NSs by 
a nonspinning BH will occur for the case that the BH 



mass is ^ 4AfQ, for the hypothetical mass and radius of 
the NS il^R ~ 1.35Mq and Rns ~ 11-12 km, respec- 
tively [a 0, S]- The tidally disrupted NSs may form a 
disk or a torus around the BH if the tidal disruption oc- 
curs outside the innermost stable circular orbit (ISGO). 
A system consisting of a rotating BH surrounded by a 
massive, hot torus has been proposed as one of the likely 
sources for the central engine of gamma-ray bursts with a 
short duration [I (hereafter SGRBs). Hence, the merger 
of a low-mass BH and its companion NS can be a candi- 
date of the central engine. According to the observational 
results by the Swift and HETE-2 satellites [13 , the to- 
tal energy of the SGRBs is larger than ^ 10''* ergs, and 
typically 10'*^-10^° ergs. The studies of hypercritical, 
neutrino-dominated accretion disks around a BH suggest 
that disk mass should be larger than ~ O.OIM© for pro- 
viding such high en ergy for generating gamma-ray via 
neutrino process [Il|,[ll|. The question is whether or not 
the mass and thermal energy of the disk (torus) are large 
enough for driving the SGRBs of the huge total energy. 
Numerical-relativity simulation plays an important role 
for answering this question. 

In the last three years, general relativistic numerical 
simulations for the inspiral and merger of BH-NS binaries 
have been performed by three groups [H, [13, [H, [l^ [l3] ■ 
However, most of these previous works should be re- 
garded as preliminary ones because the simulations were 
performed only for a short time scale: In the early works 
[T^. [l5| . the inspiral motion is followed for at most two 
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orbits. With such short-term simulations, the orbit is not 
guaranteed to be quasicircular at the onset of the merger, 
and hence, the obtained results are not likely to be very 
realistic because of the presence of the unrealistic eccen- 
tricity and incorrect approaching velocity. Moreover, by 
such incorrect conditions at the onset of the merger, the 
final outcome such as mass and spin of the BH and the 
physical condition for the disk surrounding the BH is not 
likely to be computed correctly, although a rough quali- 
tative feature of the merger mechanism and gravitational 
waveforms have been found from these works. In the lat- 
est works [l^ , the inspiral motion is followed for ~ 4 
orbits, but in each of these works, the simulation was 
performed only for one model. More systematic study 
is obviously necessary to clarify the quantitative details 
about the merger process, the final outcome, and gravi- 
tational waveforms. 

In this paper, we report our latest work in which a 
long-term simulation of BH-NS binaries is performed for 
a wide variety of BH masses, NS masses, and NS radii 
for the first time. In the present simulation, the inspi- 
ral motion is followed for 4-7 orbits. With this setting, 
the eccentricity of the last inspiral motion appears to be 
negligible and the approaching velocity at the onset of 
the merger is correctly taken into account. Furthermore, 
we systematically choose the BH and NS masses and NS 
radii for a wide range, and as a result, it becomes possi- 
ble to clarify the dependence of the merger process, final 
outcome, and gravitational waveforms on the mass ratio 
of the binary and the compactness of the NS. One draw- 
back in the present work is that the NSs are modeled by 
a simple equation of state (EOS). However, the various 
features found in this paper will qualitatively hold irre- 
spective of the EOS and thus the present work will be an 
important first step toward more detailed simulation in 
the near future in which the NSs are modeled by more 
realistic EOSs. 

The paper is organized as follows. Section II summa- 
rizes the initial conditions chosen in this paper. Section 
III briefiy describes the formulation and methods for the 
numerical simulation. Section IV presents the numeri- 
cal results of the simulation focusing primarily on the 
dependence of the merger mechanism and gravitational 
waveforms on the mass ratio and radius of NSs. Section 
V is devoted to a summary. Throughout this paper, the 
geometrical units of c = G = 1, where c and G are the 
speed of light and gravitational constant, are used, oth- 
erwise stated. The irreducible mass of a BH, rest mass of 
an NS, gravitational mass and circumferential radius of 
an NS in isolation, Arnowitt-Deser-Misner (ADM) mass 
of system, and sum of the BH and NS mass (often re- 
ferred to as the total mass) are denoted by A/bh, -W*, 
AfNS, ^nSj -M, and mo(= AfBH+AfNs), respectively. The 
mass ratio of the binary and the compactness of the NS 
are defined by Q = Mbh/Mns and C = GMns/c^^^ns, 
respectively. Note that Mbh is equal to the ADM mass 
of the BH in isolation for a nonspinning BH. Greek in- 
dices and Latin indices denote the spacetime and space 



components, respectively; Cartesian coordinates are used 
for the spatial coordinates. 



II. INITIAL CONDITION 

BH-NS binaries in a quasiequilibrium state are em- 
ployed as an initial condition of the numerical simulation. 
Following our previous papers [H, [l3| , the quasiequilib- 
rium state is computed in the moving-puncture frame- 
work p^ . [Toj . The formulation in this framework is 
slightly different from that in the excision framework 
which is adopted in most of the previous works 0, 0, 
S [2^ [2l| , although both formulations are based on the 
conformal-flatness formalism for the three-metric [l^l • In 
the present work, we adopt the same basic equations as 
those in Ref. [ij]. However, we change the condition 
for defining the center of mass of system to improve the 
quality of the quasiequilibrium (specifically, to reduce the 
orbital eccentricity). To clarify which part is changed, 
we summarize the basic equations and method for deter- 
mining the quasiequilibrium state again in the following. 
Detailed numerical solutions and their properties are pre- 
sented in an accompanied paper [23| , to which the reader 
may refer for details. 



A. Formulation 

If the orbital separation of a BH-NS binary is large 
enough, the time scale of gravitational-wave emission 
(tgw) is much longer than the orbital period (Pq). In 
the present work, we follow the inspiral motion of the 
BH-NS binaries for more than 4 orbits (typically ^ 5 or- 
bits). This implies that the initial binary separation is 
always large enough to satisfy the condition, tqw ^ ^o- 
Thus, the binaries initially given should be in a quasi- 
circular orbit (i.e., the BH and NS are approximately in 
an equilibrium in the comoving frame with an angular 
velocity il). To obtain such a state, we may assume the 
presence of a helical Killing vector around the center of 
mass of the system as follows. 



(1) 



where VL denotes the orbital angular velocity. 

We assume that the NSs have the irrotational velocity 
field because it is believed to be realistic for the BH-NS 
binaries in close orbits [l^]. For the fiuid of the irrota- 
tional velocity field, the Euler and continuity equations 
reduce to a first integral of motion and an elliptic-type 
equation for the velocity potential in the presence of the 
helical Killing vector |25|. As a result, the density profile 
and velocity field are determined by solving these hydro- 
static equations. 

For computing a solution of the geometric variables 
of a quasiequilibrium state, we employ the so-called 
conformal-fiatness formalism for the three-geometry [l^l ■ 
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In this formalism, a solution is obtained by solving 
Hamiltonian and momentum constraint equations, and 
an equation for the lapse function (a) which is derived by 
imposing the maximal slicing condition as = = dtK 
where K is the trace part of the extrinsic curvature 
Kij. These equations lead to the equations for the con- 
formal factor -0, a rescaled tracefree extrinsic curvature 
A^-' = il)^K^\ and a weighted lapse $ = aij) as 



A0 = -2^phV'' - \A,'A;r 



A ^ 



A$ = 27r<I> 



(2) 
(3) 
(4) 

where A denotes the flat Laplacian, pn = ph{au*)'^ — P, 
Ji = phau^Ui, and S = ph[{au*^)^ ~ 1] + 3P. p is the 
rest-mass density, h is the specific enthalpy defined by 
1 + e + P/ p, e is the specific internal energy, P is the 
pressure, and is the four-velocity. 

For the relation among p, s, and P, we adopt the poly- 
tropic EOS as, 



P 



Kp 



(5) 



where k is the adiabatic constant and T the adiabatic 
index for which we choose 2 in this paper. In this EOS, 
s is given by P/ (F — l)p. 

We solve the elliptic equations ©-([l]) in the moving- 
puncture framework [l^, [l^, |2^. Assuming that the 
puncture is located at rp(= Xp), we set ip and <& as 



Mp 
2rBH 



and $ = 1 - 



M<j. 



(6) 



where Mp and M$ are positive constants of mass dimen- 



sion, and tbh 



BH 



Xp). Then, substi- 



tuting Eq. ([6]) into Eqs. ^ and ([4]), elliptic equations 
for new non-singular functions (f> and 77 are derived. 

The mass parameter Mp may be arbitrarily given, and 
thus, it is appropriately chosen to obtain a desired BH 
mass. For a given value of Mp, A/$ is determined by the 
virial relation, which should hold for stationary space- 
times (e.g., Ref. [13)) as 



a$d5* = - 



^^^JdS' = 27rMo 



(7) 



where Mq is the initial ADM mass of the system. 

For a solution obtained in this formalism, a always be- 
comes negative near the puncture (approximately, inside 
apparent horizon) , but such lapse is not favorable for nu- 
merical simulation. Following previous papers [isl . 
thus, the initial condition for a is appropriately modified 
so as to satisfy the condition of a > everywhere. 

Equation ^ is rewritten by setting 



Aij (- 



-d.,S'''Wk,i+Kf^, (8) 



where Wi{= W^) denotes an auxilary three-dimensional 
function and Kf^ denotes a weighted extrinsic curvature 
associated with the linear momentum of a BH; 



K, 



^'bh 



HiPj + rijPi + {niUj - Sij)Pkn'' 



(9) 



Here, rik — — x^^/rp,ii and Pi{— P*) denotes the lin- 
ear momentum of the BH, determined from the condition 
that the total linear momentum of system is zero as 



P^ = - J^rd^X 



(10) 



The right-hand side of Eq. (|T0|) denotes the linear mo- 
mentum of the companion NS. Then, the total angular 
momentum of the system is derived from 



J= / J^ip^d^x + e^jkX^P'', 



(11) 



where we assume that the z axis is the axis of orbital 
rotation. 

Substituting Eq. ([8|) into Eq. ([3]), an elliptic equation 
for Wj is derived in the form 



(12) 



Denoting Wi — 7Bi — (x,i + Bk^ix'^) where x and Bi 
are auxiliary functions, we decompose Eq. (jl2[) into two 
linear elliptic equations 



ABi = nJi^p and Ax = —TrJix''ip' 



(13) 



To determine a BH-NS binary in quasiequilibrium 
states, in addition, we have to determine the shift vec- 
tor, /3*. The reason for this is that /3* appears in the 
hydrostatic equations [23| . In the conformal-flatness for- 
malism, the relation between and Aij is written as 

Sjkd.p" + S.kdjP" - ^S.jdkP" = . (14) 

Operating S^^di to this equation, an elliptic equation is 
derived 

A/3'' + \s'''dkdjl3^ = 2dj{ai;'^)A'^ + IG-KaJjS'^ . (15) 

For Aij in the right-hand side of this equation, we sub- 
stitute the relation of Eq. (HI [not Eq. (fT4|) ]. As a result, 
no singular term appears in the right-hand side of Eq. 
(lisp , and thus, is solved in the same manner as that 
for W^. 

In this formulation, the elliptic equations for the 
gravitational- field components, </>, ij, Bi, and /?', and 
the velocity potential have to be solved. For a numerical 
solution of them, we use the LORENE hbrary by 
which a high-precision numerical solution can be com- 
puted using the spectral method. We note that in the 
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puncture framework, we do not basically have to im- 
pose an inner boundary condition around the BH. In 
this point, the moving-puncture framework differs from 
the excision framework, and this may be a demerit of 
the moving-puncture framework because a physical con- 
dition may not be imposed for the BH. However, this 
may be also a merit of this framework because there re- 
mains a degree of freedom which can be used to adjust 
the property of the quasiequilibrium to a desired state, 
as mentioned in the following. 

The final remaining task in the moving-puncture 
framework is to determine the center of mass of sys- 
tem. The issue in this framework is that we do not have 
any natural physical condition for determining it. (By 
contrast, the condition is automatically derived in the 
excision framework Q, although it is not clear whether 
the condition is really physical and whether the result- 
ing quasiequilibrium is a quasicircular state; see, e.g., 
Ref. ^] which assesses the circularity of the quasiequi- 
librium.) In our first paper [l3| , we employed a condition 
that the dipole part of ijj at spatial infinity is zero. How- 
ever, in this method, the angular momentum derived for 
a close orbit of mo^ ^ 0.03 is by ^ 2% smaller than that 
derived by the third post-Newtonian (3PN) approxima- 
tion (Note mo = A/bh + Mns-) Because the 3PN 
approximation should be an excellent approximation of 
general relativity for describing a binary of a fairly dis- 
tant orbit as toqO « 0.03, we should consider that the 
obtained initial data deviates from the true quasicircu- 
lar state, and thus, the initial orbit would be eccentric. 
Such initial condition is not suitable for quantitatively 
accurate numerical simulation of the inspiraling BH-NS 
binaries. 

In the subsequent work 14|, we adopted a condition 
that the azimuthal component of the shift vector /J"^ at 
the location of the puncture (r = rp) is equal to —fl; i.e., 
we imposed a corotating gauge condition at the location 
of the puncture. In the following, we refer to this condi- 
tion as "(3"^ condition" . This is slightly better than the 
original condition, but the angular momentum derived 
for a close orbit of mo^l > 0.03 is still by ^ 1% smaller 
than that derived by the 3PN relation for a large mass 
ratio Q >2 (see Ref. [23^ for detailed numerical results). 
The disagreement is larger for the larger mass ratio. As 
a result of this, the initial condition is likely to deviate 
from the true quasicircular state and hence the initial 
orbital eccentricity is not also negligible (see Sec. IV C 
of Ref. [l^ for numerical evolution of such initial data) , 
in particular, for binaries of large mass ratios. This also 
suggests that the /S*^ condition is not suitable for deriving 
a realistic quasicircular state. 

If the simulations are started with a sufficiently large 
orbital separation, the eccentricity, which is initially 
~ 0.1, will decrease to ^ 0.01 within several orbits, be- 
cause gravitational radiation reaction has a strong effect 
to reduce the orbital eccentricity [lH. However, the sep- 
aration has to be large enough for the large initial ec- 
centricity to be reduced. This implies that a long-term 



simulation, which is not computationally favored, is re- 
quired. (As we illustrate in this paper, the eccentricity is 
not sufficiently suppressed to < 0.01 in ~ 5 orbits if we 
adopt the initial condition obtained in the f3^ condition.) 

In this paper, we employ a new condition in which 
the center of mass is determined in a phenomenological 
manner: We impose the condition that the total angu- 
lar momentum of the system for a given value of mo^l 
agrees with that derived by the 3PN approximation (see 
Ref. ^S^l for a similar concept). This can be achieved 
by appropriately choosing a hypothetical position of the 
center of mass (it may not be better to refer to this posi- 
tion as "the center of mass" but simply as "the rotation 
axis" ) . With this method, the total energy of the system 
does not agree completely with that derived by the 3PN 
approximation [23] . and thus, the eccentricity does not 
become zero. However, the resulting eccentricity in this 
condition is much smaller than that in the (3'^ condition 
(see Fig. [IJ, and thus, for a moderately long-term simu- 
lation (in ~ 2-3 orbits), the effect of the eccentricity is 
suppressed to an acceptable level for a scientific discus- 
sion, as shown in Sec. IIVI We refer to this condition as 
"3PN-J condition" in the following. 

Figure [T] plots evolution of the orbital separation for 
C — 0.145 and for Q = 2 and 5. Here, the separation 
is defined as the coordinate separation between the po- 
sitions of the puncture and of the maximum density of 
the NS, Tsop = la^sopl; see Eq. (PT|) . For the results with 
the initial condition derived in the P'f condition, the sep- 
aration badly oscillates with time and the amplitude of 
this oscillation is still conspicuous even after 3-4 orbital 
motions. Consequently, the eccentricity is not negligible 
even at the last orbit just before the merger. By con- 
trast, for the case that the initial condition is derived in 
the 3PN-J condition, the amplitude of the oscillation is 
much smaller and, furthermore, the amplitude of this os- 
cillation does not become conspicuous within about two 
orbital motions. Although the eccentricity does not be- 
come exactly zero, it is within an acceptable level at the 
last orbit just before the merger (see, e.g. gravitational 
waveforms shown in Sec. IIV Cp . 

We note that the coordinate separation is a gauge- 
dependent quantity and hence the discussion here is not 
based on the very physical quantity. The physical quan- 
tity such as the orbital eccentricity is not rigously ex- 
tracted from it. However, the similar quantitative fea- 
ture is seen if we plot evolution of the gravitational-wave 
frequency as a function of time (which is the physical 
quantity); the oscillation amplitude of the orbital sep- 
aration shows the magnitude of the orbital eccentricity 
at least approximately (see Sec. IIV CI for more physical 
analysis). 



B. Chosen models 

In the polytropic EOS, the adiabatic constant, k, is a 
free parameter, and thus, physical units such as mass. 
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(a) (b) 

FIG. 1: Evolution of orbital separation (a) for model M20.145 and (b) for model M50.145. The solid and dotted curves denote 
the results with the initial conditions obtained in the 3PN-J and (5^ conditions, respectively. The dashed curve in panel (a) 
denote the result for model M20.145N. To align the curves at the onset of the merger, the time is appropriately shifted for 
the results of M20.145N, M20.145b, and M50.145b. Note that the merger sets in when the orbital separation becomes ~ SMq. 
The binaries for models M20.145 and M50.145 spend in the inspiral phase for ~ 5 and 7.5 orbits, respectively (see Sec. IV), 
whereas those for models M20.145b and M50.145b spend only for ~ 4 and 5.25 orbits, respectively. 



TABLE I: Key parameters for the initial conditions adopted for the numerical simulation in units of k = 1. The mass ratio 
(Q = Mbh/Mns), BH mass (Mbh), mass parameter of the puncture {Mp), rest mass of the NS (M*), mass (Mns) and 
compactness (C = Mns / ^?ns : ^ns is the circumferential radius) of the NS when it is in isolation, maximum density of the NS 
(pmax), ADM mass of the system (Mo), total angular momentum (Jo) in units of Mq , orbital period (Pq) in units of Mo, and 
orbital angular velocity (f2o) in units of utIq^ . The first and second numerical values described for the model name (the first 
column) denote the values of Q and C, respectively. The first 10 models are computed in the 3PN-J condition and the last 4 
are in the /S*^ condition. 



3PN-J 
Model 


Q 


A/bh 


A'/p 


M, 


A/ns 


A/ns / -Rns 


Pmax 


Mo 


Jo /MS 


Pq/A/q 


mo^o 


M15.145 


1.5 


0.2093 


0.2043 


0.1500 


0.1395 


0.145 


0.1262 


0.3452 


0.8393 


212.5 


0.02987 


M20.145 


2 


0.2790 


0.2736 


0.1500 


0.1395 


0.145 


0.1261 


0.4146 


0.8603 


213.7 


0.02968 


M20.145N 


2 


0.2790 


0.2732 


0.1500 


0.1395 


0.145 


0.1260 


0.4143 


0.8418 


191.8 


0.03309 


M20.160 


2 


0.2957 


0.2900 


0.1600 


0.1478 


0.160 


0.1512 


0.4394 


0.8608 


214.4 


0.02958 


M20.178 


2 


0.3119 


0.3059 


0.1700 


0.1560 


0.178 


0.1890 


0.4635 


0.8582 


211.4 


0.03001 


M30.145 


3 


0.4185 


0.4121 


0.1500 


0.1395 


0.145 


0.1255 


0.5534 


0.7091 


192.2 


0.03298 


M30.160 


3 


0.4435 


0.4367 


0.1600 


0.1478 


0.160 


0.1504 


0.5864 


0.7096 


192.9 


0.03285 


M30.178 


3 


0.4679 


0.4608 


0.1700 


0.1560 


0.178 


0.1878 


0.6187 


0.7103 


193.8 


0.03269 


M40.145 


4 


0.5580 


0.5512 


0.1500 


0.1395 


0.145 


0.1252 


0.6926 


0.6042 


192.1 


0.03294 


M50.145 


5 


0.6975 


0.6905 


0.1500 


0.1395 


0.145 


0.1249 


0.8318 


0.5238 


191.8 


0.03296 


Model 


Q 


Mbh 


Mp 


M* 


Mns 


Mns//?ns 


Pmax 


Mo 


Jo /MS 


Po/Mo 


mo^o 


M20.145b 


2 


0.2790 


0.2737 


0.1500 


0.1395 


0.145 


0.1261 


0.4144 


0.8462 


211.0 


0.03007 


M30.145b 


3 


0.4185 


0.4121 


0.1500 


0.1395 


0.145 


0.1256 


0.5531 


0.6960 


189.5 


0.03345 


M40.145b 


4 


0.5580 


0.5512 


0.1500 


0.1395 


0.145 


0.1252 


0.6922 


0.5905 


188.9 


0.03351 


M50.145b 


5 


0.6975 


0.6905 


0.1500 


0.1395 


0.145 


0.1250 


0.8315 


0.5134 


189.1 


0.03345 



radius, and time can be rescaled arbitrarily by simply 
changing the value if k: i.e., when a numerical result for 
a particular value (say k = ki) is obtained, we can also 
obtain the numerical results of these quantities for n — k,2 
by simply rescaling them by a factor of (k2/ki)^^^''^~^^. 
This implies that k can be completely scaled out of the 
problem. In this paper, we present the results in units 



of K = 1 (and c ^ G — 1}, because such units are popu- 
lar in other groups [ll,[l3|- the polytropic EOS, the 
nondimensional quantities such as C = GA/ns/c^-RnSi 
Mofl, MoK-i/2(r-i), i?NsK-^/2(^~^\ and mass ratio (Q) 
are unchanged irrespective of the value of k and have an 
invariant meaning. 

In the present work, the numerical simulation is per- 
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0.13 0.14 0.15 0.16 0.17 0.18 0.19 

FIG. 2: Initial conditions in the parameter space of (C, Q) are 
plotted by the solid circles. The meaning of the dashed curve 
is as follows [3: For the binaries located above the dashed 
curve, mass-shedding (MS) does not occur until the binaries 
reach the ISCO, whereas for the binaries below the dashed 
curve, the mass-shedding will occur for the NS due to the 
tidal force of the BH before the ISCO is reached. 



formed for a wide variety of initial conditions (see Ta- 
ble I), but for restricting that the BH spin is zero. 
We characterize the BH-NS binaries by the mass ra- 
tio, Q = A/bh/A/nSj and the compactness of the NS, 
C = M^s/Rm- (Note that in Refs. [13, [ll], we use 
q = 1/Q instead of Q to specify the model.) The mass 
ratio is chosen in the range between 1.5 and 5.0, and the 
compactness of the NS is in the range, 0.145-0.178. We 
note that the typical mass of the NS in nature is 1.3- 
1.4Mns (sS*! and the likely lower bound of the BH mass 
is ~ 2Mq. This implies that Q should be chosen to be 
larger than ^ 1.5. The circumferential radius of an NS 
of Mns = 1.35M0 and C = 0.145, 0.160, and 0.178 is 
13.8, 12.5, and 11.2 km, which are reasonable values for 
modeling the NSs. 

Figure [2] shows the initial conditions in the parameter 
space of (C, Q). The meaning of the dashed curve, de- 
rived in Ref. is as follows: For the binaries located 
above the dashed curve, mass-shedding of the NS by the 
tidal effect of the companion BH does not occur until the 
binary reaches the ISCO, whereas for the binary shown 
below the dashed curve, the mass-shedding will occur for 
the NS before the ISCO is reached and in such a system, 
tidal disruption may subsequently occur. Thus, many of 
the NSs in the chosen binary system will be subject to 
the tidal effect of the companion BH in a close orbit, but 
for some of them (e.g., {C,Q) = (0.145,5) and (0.178, 
3)), the tidal effect is unlikely to play an essential role in 
the inspiral phase. We also note that an NS which is in 
the mass-shedding is not always tidally disrupted imme- 
diately, although the mass-shedding NS is the candidate 
to be tidally disrupted. Namely, the condition for induc- 
ing the tidal disruption is in general more restricted than 
that for inducing the mass-shedding (see Sec. IV). 

In the present work, we prepare quasiequilibrium states 
with moflo « 0.030 for Q = 1.5 and 2, and mo^lo ~ 0.033 



for Q = 2-5, where flo is the orbital angular velocity. We 
choose these values of JIq so as for the BH-NS binaries 
to experience more than 4 orbits before the onset of the 
merger. In such a long-term evolution, the eccentricity 
which presents at t = decreases to an acceptable level 
during the inspiral phase, and also the nonzero radial 
velocity associated with the gravitational radiation reac- 
tion is approximately taken into account from the late 
inspiral phase. 

In Table I, several key quantities for the quasiequilib- 
rium states adopted in this paper are listed. Specifically, 
we prepare 10 models. The first two and last three nu- 
merical numbers described for the model name denote 
the values of Q and C, respectively. For comparison, the 
same quantities are also listed for 4 selected quasiequilib- 
rium states obtained in the Z?*^ condition. We find that 
the angular momentum of them is by ~ 2% smaller than 
that of the quasiequilibrium states obtained in the 3PN- J 
condition for moflo « 0.033. 

III. PREPARATION FOR NUMERICAL 
SIMULATION 

A. Brief summary of formulation and methods 

The numerical simulations are performed using a code 
SACRA recently developed in our group [l^. The details 
of the chosen scheme, formulation, gauge condition, and 
methods for the analysis are described in Ref. [l6| to 
which the reader may refer. Reference [lB| also shows 
that SACRA can successfully simulate the inspiral and 
merger of BH-BH, NS-NS, and BH-NS binaries. 

In SACRA, the Einstein equations are solved in a 
moving-puncture version [l^, [3, HI] of the Baumgarte- 
Shapiro-Shibata-Nakamura formalism [s^]. Specifically, 
we evolve W = 7"^/^ [3^, conformal three- metric, 
7y = 7~^^'^7ij, the tracefree extrinsic curvature, Aij = 
^~^/^{Kij — K^ij/3), the trace part of Kij, K, and a 
three auxiliary variable, F* = —djj^^ (or Fi = djjij). 
Here, 7^ is the three-metric, Kij the extrinsic curvature, 
and 7 = det(7y ). In the numerical simulation, a fourth- 
order finite differencing scheme in space and time is used 
implementing an adaptive mesh refinement (AMR) al- 
gorithm (at refinement boundaries, a second-order in- 
terpolation scheme is partly adopted). The advection 
term such as P^dijjk is evaluated by a fourth-order non- 
centered finite-difference, as proposed in Ref. [26j . 

The moving-puncture approach is adopted for evolving 
the BH [T^ (see also Ref. [1^); i.e., we adopt one of the 
moving-puncture gauge conditions as (26j 

{dt- f3^d.j)p' = 0.75B\ (16) 
{dt - f3'd,)B' - (dt - f3'd,)r - ij,B\ (17) 

where i?* is an auxiliary variable and rjs is an arbitrary 
constant. In the present paper, we set rjs — 0.414 in units 
ofK = c = G = l. 
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TABLE II: Parameters of the grid structure for the numerical simulation with our AMR algorithm. In the column named 
"Levels", the number of total refinement levels is shown (in the bracket, the numbers of coarser and finer levels are specified; 
see Ref. ^ for the definition of the coarser and finer levels). Aa;( — /17) is the minimum grid spacing for — 36, -Rdiam 
the coordinate length of semi-major diameter of the NS, L the location of outer boundaries along each axis, Ao(= n/Qo) the 
gravitational wavelength at t = 0, and Axgw the grid spacing, by which gravitational waves are extracted for A'^ = 36. Note 
that Mbh denotes the BH mass (the irreducible mass) at t = 0. The grid structures for models M20.145b-M50.145b are the 
same as those for models M20.145-M50.145, respectively. For the simulations with A' 7^ 36, the size of each domain (and hence 
L) is unchanged, and thus, the grid spacing hi is simply changed. 



Run 


Levels 


Ax /Mo {Ax /Mbh) 


i?diam/Ax 


L/Mo (L/Ao) 




JMo 


M15.145 


8 (4+4) 


0.0407 (0.0672) 


112 


187.7 (1.77) 


1.30- 


-5.21 


M20.145 


8 (4+4) 


0.0377 (0.0560) 


99.3 


173.7 (1.63) 


1.21- 


-4.82 


M20.145N 


8 (4+4) 


0.0377 (0.0560) 


99.4 


173.8 (1.81) 


1.21- 


-4.83 


M20.160 


8 (4+4) 


0.0356 (0.0528) 


92.9 


163.9 (1.53) 


1.14- 


-4.55 


M20.178 


8 (4+4) 


0.0324 (0.0481) 


89.1 


149.1 (1.41) 


1.04- 


-4.14 


M30.145 


8 (4+4) 


0.0305 (0.0403) 


90.0 


140.5 (1.46) 


0.98- 


-3.90 


M30.160 


8 (3+5) 


0.0266 (0.0352) 


91.2 


122.8 (1.27) 


0.85- 


-3.41 


M30.178 


8 (3+5) 


0.0253 (0.0334) 


84.1 


116.3 (1.20) 


0.81- 


-3.23 


M40.145 


8 (3+5) 


0.0244 (0.0302) 


89.0 


112.3 (1.17) 


0.78- 


-3.12 


M50.145 


8 (3+5) 


0.0203 (0.0242) 


88.4 


93.5 (0.98) 


0.65- 


-2.60 



For the numerical hydrodynamics, we evolve p* = 
pait'e^*^, Ui = hui, and e* = phau* — P/{pau^). To han- 
dle the advection terms in the hydrodynamic equations, 
a high-resolution central scheme [s^l is adopted with a 
third-order piecewise parabolic interpolation and with a 
steep min-mod limiter in which the limiter parameter b 
is set to be 3 (see appendix A of Ref. [111). We adopt 
the F-law EOS in the simulation as 

P = (F - l)p£, (18) 

where F = 2. 

Properties of the BHs such as mass and spin are de- 
termined by analyzing area and circumferential radii of 
apparent horizons. A numerical scheme of our apparent 
horizon finder is described in Refs. (TgI. [39j. 

Gravitational waves are computed by extracting the 
outgoing part of the Newman-Penrose quantity (the so- 
called ^"4). The extraction of is carried out for several 
constant coordinate radii, r m SO-IOOMq. The plus and 
cross modes of gravitational waves are obtained by per- 
forming time integration of ^'4 twice, with appropriate 
choice of integration constants and subtraction of un- 
physical drift which is caused primarily by the drift of 
the center of mass of the system. (Because we extract 5*4 
for fixed, finite coordinate radii, the drift of the center of 
mass spuriously affects gravitational waveforms.) Specif- 
ically, whenever the time integration is performed, we 
subtract a function of the form a2t^ + ait + ao where a^- 
a2 denote constants which are determined by the least- 
square fitting of the numerical data. 

We compute the modes of 2 < Z < 4 and \m\ < I for '^4, 
and found that the quadrupole mode of (Z, |m|) = (2,2) 
is always dominant, but (Z,|m|) = (3,3), (4, 4), and (2, 
1) modes also contribute to the energy and angular mo- 
mentum dissipation by more than 1% for some of models 



(in particular for large values of Q; cf. Table IV). 

We also estimate the kick velocity from the linear mo- 
mentum fiux of gravitational waves. The linear momen- 
tum flux dPi /dt is cornputed by the same method as that 
given in, e.g., Refs. [l^, H^l- Specifically, the couphng 
terms between (Z, m) = (2, ±2) and (3, +3) modes, be- 
tween (2, ±2) and (2, +1) modes, and between (3, ±3) 
and (4, +4) modes contribute primarily to the linear mo- 
mentum fiux [4]| . From the total linear momentum dis- 
sipated by gravitational waves, 

AP../f<«, (19) 

the kick velocity is defined by APi/Mo where Mo is the 
initial ADM mass of the system. 

B. Setting grids for AMR scheme 

The Einstein and hydrodynamic equations are solved 
in an AMR algorithm described in Ref. [l^. In the 
present work, we prepare 8 refinement-level domains of 
different grid resolutions and domain sizes. Each domain 
is composed of the uniform vertex-center-grid with grid 
number {2N + 1, 27V + 1, iV + 1) for {x, y, z) where N is 
a constant and is always chosen to be 24, 30, and 36 to 
check convergence of numerical results. The equatorial 
plane symmetry with respect to the z = plane is as- 
sumed. The length of a side for the largest domain is 
denoted by 2L, and the grid spacing for each domain is 
hi = L/N/2^ for I = 0-7. As described in Ref. [16], the 
regions around a BH and an NS are covered by "finer" 
domains which move with the BH and NS. On the other 
hand, "coarser" domains which cover wider regions do 
not move and their center is fixed approximately to be 
the center of mass of the system. 
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Table II lists the parameters for the grid structure in 
our AMR scheme. Aq denotes the wavelength of gravi- 
tational waves at t — (Aq = tt/Hq). For all the cases, 
L is chosen to be l-2Ao, implying that the outer bound- 
aries are located in a wave zone. The NS is covered by 
the finest and second-finest domains and the coordinate 
radius of the apparent horizon for the BH is always cov- 
ered by more than 15 grid points for iV = 36. Because 
the numerical results (such as the time spent in the in- 
spiral phase, the mass and spin for the final state of the 
BH, and total energy radiated by gravitational waves) for 
iV > 30 depend only weakly on the grid resolution, we 
conclude that the convergence is approximately achieved 
for = 36 (except for the rest mass of disks which are 
formed for the model with small values of C and Q; see 
Sec. HVH) . 

For = 36, the total memory required for the simu- 
lation with 13 domains is about 5 GBytes. We perform 
all the simulations using personal computers of 8 GBytes 
memory and 2-8 processors (in one job, only two pro- 
cessors are used with an OpenMP library). The typical 
computation time for one model with A^ = 36 is 5-7 
weeks on the personal computers of clock speed 3 GHz. 



C. Atmosphere 

Because any conservation scheme of hydrodynamics is 
unable to evolve a vacuum, we have to introduce an ar- 
tificial atmosphere outside the NSs. The density of the 
atmosphere has to be small enough to exclude its spuri- 
ous effect to the orbital motion of the BH and NS, and to 
avoid overestimation of the total rest mass for a disk sur- 
rounding the BH, which may be formed after the merger 
if the NS is tidal disrupted. We initially assign a small 
rest-mass density for the atmosphere as follows 



Patmo 



Pcrit 
Pcrit 



r < To 

gl-r/ro > 



(20) 



where Tq is a constant chosen to be 5L/32. We choose 
Pcrit = Pmax X 10~^ whcrc Pmax is the maximum rest- 
mass density of the NS initially given. With such choice, 
the total amount of the rest mass in the atmosphere is 
less than 10~^M*. Accretion of the atmosphere of such 
small mass onto the BH and NS plays a negligible role for 
their orbital evolution. As we state in Sec. I, one of the 
important tasks for the merger simulation of the BH-NS 
binaries is to determine the rest mass of disks surround- 
ing a BH formed after the merger. In the simulation, we 
pay attention only to the case that the disk mass is larger 
than IQ-'^Af*. 

During the evolution, we also adopt an artificial treat- 
ment for the low-density region in the following manner: 
(i) If the density is smaller than Patmo, we set p — Patmo 
and Ui = 0. Then, the specific internal energy s is set to 
be Kp^^^/{r — 1). (ii) Even if the density is larger than 
Patmo, we reduce the specific momentum hui by a factor 



of 1 — exp[— p/pcrit]; i-e., for a fluid of density ^ 5pcrit, 
the specific momentum is artificially reduced. The rea- 
son that we adopt these treatments is that a numerical 
instability resulting in a negative density or pressure of- 
ten happens accidentally for the low-density region in the 
absence of the artificial treatment. However, by limiting 
the unphysical growth of the specific momentum, as de- 
scribed above, such instabilities are excluded. For some 
models, we checked whether the magnitude of pcnt af- 
fects the numerical result, but as long as the small value 
is chosen for it as in the present work, dependence of the 
numerical results on the magnitude of pcrit is quite tiny. 



IV. NUMERICAL RESULTS 

A. Orbital evolution, tidal disruption, and disk 
mass 

Figure [3] plots evolution of the coordinate separation 
between the BH and NS, a;L,„ for models M20.145, 



M40.145, and M50.145. 



^scp 



scp ' 



scp 



is defined by 



(21) 



where xj^g and denote the positions of the maximum 
rest-mass density of the NS and of the puncture, respec- 
tively. This figure illustrates that the binaries are in a 
slightly eccentric orbit for the first ^ 2 orbits because 
the strictly circular orbit is not provided initially. Also, 
in the first ~ 2 orbits, decrease rate of the orbital sep- 
aration due to gravitational radiation reaction is not as 
large as that predicted by the PN theory, in particular 
for larger values of Q. However, because the eccentric- 
ity decreases due to the gravitational radiation reaction 
during the evolution and also the initial eccentricity is 
not very large (see, e.g.. Fig. [1]), the orbit approaches 
approximately to a quasicircular orbit after a few orbits. 
The resulting final 2-3 orbits before the onset of merger 
appear to be close to a quasicircular one. This behavior 
is much better than that obtained when an initial condi- 
tion computed in the f]"^ condition is adopted (see, e.g. 
Fig. 15 of Ref. 13). 

For models M20.145, M40.145, and M50.145, the bi- 
naries spend in the inspiral phase for ~ 4.8, 6.5, and 7.8 
orbits, respectively. For the larger values of Q with an ap- 
proximately fixed value of mo^o, the number of the inspi- 
ral orbit is larger, because the luminosity of gravitational 
waves is approximately proportional to Q^/ {1 + QY, and 
as a result, the binary evolves more slowly for the larger 
values of Q. 

Figures [4]-[6] plot late-time evolution of the rest-mass 
density contour curves and density contrasts as well as 
the location of the BHs for models M20.145, M40.145, 
and M50.145, respectively. For model M20.145, the NS 
is tidally disrupted in the late inspiral phase (see the 
first panel of Fig. Subsequently the material of the 
NS forms a one-armed spiral arm around its compan- 
ion BH (second panel of Fig. 2]). In the spiral arm, a 
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FIG. 4: Snapshots of the density contour curves and density contrasts as weU as the location of the BH in the merger and 
ringdown phases for model M20.145. The contour curves are plotted for p — 10~* where i = 2,3,4,5 in the first four panels, 
whereas in the last two, p = 10~' where i = 3, 4, and 5. The first panel denotes the state just at the onset of the merger. The 
filled circles show the region inside the apparent horizon. Note that the maximum value of p for the NS in the inspiral phase 
is ?a 0.126. 



transport process of angular momentum from its inner 
to outer region is likely to work efficiently. Then, the 
spiral arm winds around the BH (third panel of Fig. [4]), 
and most of the fluid elements, which do not have an- 
gular momentum large enough to have an orbit around 
the BH, fall into the BH. In the first - 200Mo after the 
tidal disruption, ^ 98% of the material falls into the BH 
(see Fig. [T]). However, a small fraction of the fluid ele- 
ments obtain angular momentum large enough to escape 
the capture by the BH, and form a disk around the BH. 
For t — tmcrgcr ^ SOOMq, whcrc ^merger approximately 
denotes the time at which the merger sets in, accretion 
rate decreases, and then, the disk relaxes to a quasisteady 
state (fourth-sixth panels). For model M20.145, the rest 
mass of the disk is ~ O.OIM* at t — tmorgor ~ lOOOAfo 
for = 36 (cf. Fig. [7]). For a hypothetical value of 
A/ns — 1.35Mo, lOOOMo is approximately equal to 20 
ms, and thus, the lifetime of the formed disk is likely 
to be much longer than 20 ms. Also, for a hypotheti- 
cal value of Mns = 1.35M0, p — 10^"' in the units of 
c = G = K = 1 corresponds to p « 6.0 x 10^^ g/cm^. 
Thus, for that hypothetical mass, the rest-mass density 
of the disk is high as ^ 10^^-10"'^^ g/cm"^. Evolution and 



the final outcome for model M15.145 are similar to those 
for model M20.145, although the disk mass is by a factor 
of ~ 2 larger due to the smaller value of Q. 

The NS for model M40.145 is also subject to tidal de- 
formation and mass shedding by the tidal effects of the 
companion BH in the late inspiral phase, as indicated in 
Fig. [Hand in the first panel of Fig. [51 However, the tidal 
effects in this binary become important only for the in- 
spiral orbits close to the ISCO. Because the approaching 
velocity is a substantial fraction of the orbital velocity at 
such close orbits, the time scale available for the NS to 
be tidally deformed before the onset of the merger is too 
short to efficiently transport angular momentum inside 
the NS. As a result, one-armed spiral arm is formed in 
a less conspicuous manner than that for model M20.145 
(see the second panel of Fig. [5]), although the NS is highly 
elongated at the merger. Rather, most of the material of 
the NS falls into the BH in a short time scale ~ 200Mo 
(see the third panel of Fig. [5]). A tiny fraction of the 
material still spreads outward during the merger phase 
(see the fourth panel of Fig. [5]), but specific angular mo- 
mentum for such material is not large enough to escape 
the capture by the BH. In this case, more than 99.999% 
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FIG. 5: The same as Fig.|4]but for model M40.145. 



FIG. 6: The same as Fig.|4]but for model M50.145. 
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of the material is eventually swallowed by the BH (cf. 
Fig. ED. 

For model M50.145, even the mass shedding does not 
occur in the inspiral phase, as indicated in Fig. (5] Dur- 
ing the merger, the NS is deformed by the tidal field 
of its companion BH (see the first and second panels of 
Fig. El) I but spiral arm is not formed nor angular momen- 
tum transport work. As a result, nearly all the materials 
are swallowed by the BH in a short time scale of ^ IOOMq 
(see the third and fourth panels of Fig. [5]) , and the final 
outcome is a rotating BH approximately in a vacuum 
spacetime. 

To clarify the infall process of the material into the 
companion BH, we plot evolution of the rest mass of the 
material located outside the apparent horizon, Mr>r ahj 
for several models in Fig. [T] Here, Mr>rAH is defined by 

Mr>rAn = / p*d^x, (22) 

•Jr>rAii 

and rAnid, ^p) denotes the radius of the apparent horizon 
for given angular coordinates. 

Figure [7(a) plots MryrAu ^ ^ function of i — inicigci 
where imorgo denotes the approximate onset time of 
the merger, for models M15.145, M20.145, M30.145, 
M40.145, and M50.145. This shows that (i) for Q < 2, a 
disk of rest mass ~ 0.01- 0.02M* is formed and the life- 
time of the disk is much longer than the dynamical time 
scale; (ii) for Q > 3, approximately all the materials are 
swallowed by the BH in ~ SOOMq (the rest mass of the 
material located outside the apparent horizon is less than 
10^^ Af* at the final stage; see also Table HI). The time 
scale for the infalling is shorter for the larger values of Q. 
These facts hold irrespective of the initial orbital sepa- 
ration and grid resolution, as long as Q > 3. Remember 
the fact that the compactness of the NSs for these mod- 
els is 0.145, which is a relatively small value; the typical 
compactness for the NS of mass 1.3-1.4Mq is between 
~ 0.14 and ~ 0.21 according to theories of high-density 
matter [42|. Because the less compact NS (i.e., the NS of 
larger radius) is more subject to disk formation, we con- 
clude that the disk (or torus) mass around a BH formed 
after the merger is negligible, for the mass ratio Q > 3. 
By contrast, if the mass ratio is smaller than ~ 2 and the 
compactness of the NS is relatively small as C = 0.145, 
a disk of mass ~ 0.01-0.02M* may be formed. (We note 
that in the presence of BH spin, this conclusion changes; 
work in progress in our group.) 

To show the further evidence that the disk is really 
formed for Q = 2 and C = 0.145, we generate Fig. [71[b). 
In this figure, we plot Mr^rAa ^ ^ function of time for 
model M20.145, for model M20.145N for which the initial 
orbital separation is smaller than that of model M20.145, 
and for model M20.145b for which the initial condition 
is computed in the (3"^ condition. This figure shows that 
the disk mass at t — tmorgor ~ IOOOMq is ~ O.OIM* ir- 
respective of the initial conditions for a given grid res- 
olution {N = 36). This fact indicates that the initial 
orbital separation is large enough to exclude spurious ef- 



fects associated with noncircularity of the initial condi- 
tion to the resulting disk mass. We note that in the 
previous early-stage works [ij, [l^ , the simulations were 
performed from the initial conditions of small orbital sep- 
arations, and consequently, it was found that the disk 
mass depends strongly on the initial separation, failing 
to draw the definitive conclusion about the disk mass. 
This drawback is overcome in this work. 

To show dependence of A/r>rAH ^ function of t on 
the compactness of the NSs, in Fig. EKc), we compare 
the results for models M20.145, M20.160, and M20.178, 
and for models M30.145, M30.160, and M30.178. As this 
figure shows, the disk mass systematically and steeply de- 
creases with the increase of the compactness C for Q = 2; 
the disk mass for model M20.160 is by a factor of ~ 15 
smaller than that for model M20.145. This dependence 
is simply caused by the fact that the tidal disruption of 
more compact NSs occurs for an orbit closer to the ISCO, 
suppressing disk formation. This result implies that even 
for a small mass ratio Q = 2, disks are not formed if the 
radius of the NS is not very large. Note that for a hy- 
pothetical mass of A/ns = 1.35A/0 for model M20.160, 
the circumferential radius of the NS is i?Ns = 12-5 km, 
which is not a large value because nuclear theories pre- 
dict it in the range ~ 10-15 km Nevertheless, the 
disk is not formed. This implies that for a typical NS 
of Mq = 1.35Af0 and i?NS = 11-12 km, the disk is 
formed only for a highly restricted case, Q < 2, i.e., 
Mbu < 2.7 M pj. This conclusion agrees with a conjec- 
ture by Miller W3| , and is also consistent with the results 
by Duez et al. |l7| in which they show that the disk mass 
is at most ^ OMM^ for a compact NS of C = 0.174 with 
the most optimistic mass ratio Q = 1. The present re- 
sult also suggests that most of BH-NS binaries may not 
be a promising candidate for the central engine of SGRBs 
11, 12], unless the radius of the NS is fairly large ^ 14 
km or the BH has a spin. 

Figure [TKc) also compares Afr>rAH a function of 
t for models M30.145, M30.160, and M30.178. For aU 
these cases, approximately all the materials of the NS 
eventually fall into the BH in ~ 500Afo after the onset 
of the merger. However, the merger process and subse- 
quent infalling process into the BH depend strongly on 
the compactness of the NSs. For model M30.145, the 
mass-shedding and subsequent tidal disruption occur be- 
fore the binary reaches the ISCO, as in model M20.145 
(see also Fig. 18 of Ref. jl^ for which essentially the 
same result as for model M30.145 is shown). As a result, 
a spiral arm is formed around the companion BH and a 
fraction of the material spreads outward. Then, a disk 
of rest mass O.IA/* surrounding the BH is transiently 
formed, although a large fraction of the material is swal- 
lowed by the BH in ~ 50A^o after the onset of the merger. 
Subsequently, the material gradually accretes from the 
disk into the BH in the time duration of ~ 15QMq. Dur- 
ing this phase, the disk mass decreases from ~ 0.1 Af* 
to 0.02Af*, and thus, for the first 150Afo after its 
formation, a massive disk is present. However, the disk 
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FIG. 7: Evolution of the rest mass of the material located outside the apparent horizon (a) for models M15.145, M20.145, 
M30.145, M40.145, and M50.145 with iV = 36, (b) for model M20.145, M20.145N, and M20.145b with N = 36, (c) for 
models M20.145, M20.160, M20.178, M30.145, M30.160, and M30.178 with N = 36, and (d) for models M15.145, M20.145, 
M20.160, and M20.178 with = 36 (solid curves) and iV = 30 (dotted curves). For a hypothetical value of Mns = 1.35M0, 
lOOAfo ~ 1.66, 2.00, 2.66, 3.33, and 3.99 ms for Q = 1.5, 2, 3, 4, and 5, respectively. 



material does not have sufficiently large specific angular 
momentum for keeping orbits around the BH, and even- 
tually falls into the BH in a runaway manner. In the end, 
more than 99.999% of the material is swallowed by the 
BH. 

For model M30.160, a disk is formed around the BH 
transiently. However, its mass is much smaller than that 
for model M30.145 because the tidal disruption occurs 
at an orbit close to the ISCO, as in model M40.145. For 
model M30.178, a disk is not formed around the BH even 
transiently, as in model M50.145. The reason is that 
the NS in this model is not tidally disrupted before the 
binary reaches the ISCO. For these two cases, more than 
99.999% of the material is swallowed into the BH within 
a short duration of 200-300Mo. 

Before closing this section, we note that the final disk 
mass depends very weakly on the grid resolution for mod- 
els M30.145, M30.160, M30.178, M20.160, and M20.178, 
whereas for models M15.145 and M20.145 for which the 
disk mass is ^ O.OIM*, the final disk mass increases with 
improving the grid resolution [see Fig. [71[d)]. This im- 
plies that (i) for the case that the massive disk is not 
formed, our conclusion is based on the convergent result. 



whereas (ii) for the case that a disk of mass ^ O.OIM* is 
formed, the results with A'' = 36 should be regarded as a 
lower-bound for the disk mass, and the disk mass would 
be larger than O.OIM* for model M20.145 and 0.02M* 
for model M15.145. This dependence on the grid res- 
olution results from the fact that with the poorer grid 
resolutions, numerical dissipation of angular momentum 
is larger, increasing an amount of the material which falls 
into the BH. However, this systematic dependence shows 
that the disk of mass > O.OIM* is indeed formed. 



B. Black hole mass and spin after merger 

Figure [8] plots Mi„./Mo, Cp/Ce, and Ce/47rMo of a BH 
as functions of time for models M20.145 and M40.145 
(similar behavior is found for other models). Here, Cp 
and Ce are polar and equatorial circumferential radii of 
the BH. An irreducible mass, Mi„, of the BH is defined 
by the area of apparent horizon, Aah, 
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FIG. 8: Evolution of Cp/Ce, Ce/47rMo, and Mirr/Mo as func- 
tions of time (a) for models M20.145 and (b) M40.145. 



Ce/47r is equal to the BH mass in stationary vacuum 
spacetimes of a BH. We follow its evolution, assuming 
that it is approximately equal to the BH mass even in 
the dynamical spacetime. 

Figure [8] shows that the values of these three quantities 
remain approximately constant before the onset of the 
merger (more specifically, before the material of the NS 
falls into the companion BH). Because the BH is not spin- 
ning initially, the hypothetical "BH mass" , Ce/47r, should 
be approximately equal to the irreducible mass Mi„, and 
Cp/Ce is approximately equal to unity. These hold ex- 
cept for small numerical error of magnitude ~ 1%. After 
the onset of tidal disruption, Ce/47r and Mi„ quickly in- 
crease as the material of the NS falls into the BH, and 
finally, they approximately reach constants. By contrast, 
Cp/ Ce decreases due to spin- up of the BH caused by mass 
accretion. Because of the presence of the BH spin, Ce/47r 
becomes unequal to the irreducible mass after the onset 
of the merger. 

In addition to Ce/47r, the mass of a BH formed after 
merger may be estimated approximately by evaluating 
the total energy dissipated by gravitational waves, AE, 
and the baryon rest mass of disks surrounding the BH 
from an approximate relation of energy conservation as 



We note that in this formula, we ignore the binding en- 
ergy between the BH and surrounding material. Thus, 
Mbh.i is likely to give a slightly overestimated value for 
the true BH mass. 

The values of AE, Mr^r ah^ -^BH.f, and Ce/47r are 
listed for all the models chosen in this paper in Table IIIII 
and IV. In Table llllj the results in the end of the simu- 
lation for = 36 are presented, whereas Table IV lists 
the results for energy and angular momentum radiated 
by gravitational waves. Here, the values of AE and A J 
depend on the radii of the wave extraction by 1-3% for 
N = 36. In addition, these systematically increase with 
improving the grid resolution. Thus, we infer these quan- 
tities by an extrapolation of the data for = 30 and 36, 
which carried out assuming the third-order convergence. 
(We note that the Einstein and hydrodynamic equations 
are solved in the fourth and third order schemes.) 

Table Hh] shows that Mbh,! agrees with Ce/47r within 
^ 0.5% error, but A/bhj is systematically larger than 
Ce/47r. The likely reason that AE, which is used in 
computing Men.f , is slightly underestimated for A^ = 36 
because of numerical dissipation of gravitational wave 
amplitude: Indeed, numerical results for the energy and 
angular momentum radiated by gravitational waves in- 
crease with improving the grid resolution as mentioned 
above, and for the extrapolated results shown in Table 
IV another conservation relations such as 



^ + Mr>rAH +AE = Mo 

47r 



(25) 



holds in a good manner within the numerical error of 
- 0.1-0.2%. 

As described in Ref. [l3| , there are at least three meth- 
ods for approximately estimating the final BH spin. In 
this paper, we use the following methods. In the first 
method, we approximately estimate the mass and spin 
of the BH by conservation laws. Namely, we determine 
them by subtracting total energy and angular momentum 
dissipated by gravitational waves, and rest mass and an- 
gular momentum of disks surrounding the BH from the 
initial ADM mass and angular momentum, respectively. 
As shown in Eq. the BH mass is estimated to give 

-^BH.f- In the same manner, angular momentum of the 
BH may be estimated by 



Jbhj — Jo — J-i 



r>rAH 



-A J, 



(26) 



where A J is total angular momentum radiated by grav- 
itational waves and Jr>rA-a is angular momentum of the 
material located outside the apparent horizon, defined by 



>''AH 



(27) 



'•>rAH 



M, 



BH,f 



Mo - M, 



r>rAii 



AE. 



(24) 



Here, u^p ~ {x ~ xp)uy + {y ~ yp)ux and {xp,yp) denote 
the position of the puncture. Jr>rAH exactly gives the an- 
gular momentum of the material in the axisymmetric and 
stationary spacetime. In the late phase of the merger, 
the spacetime relaxes to a quasistationary and nearly 
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axisymmetric state. Thus, we may expect that Jr>rAH 
will provide an approximate magnitude of the angular 
momentum of disks. From Jbh,! and Mbu,!, we define 
a nondimensional spin parameter by an = JBH,f/-^^BH f 
(see Table III). 

In the second method, the spin is determined from the 
following geometric quantities of the apparent horizon; 
Ce/47r and Mi„. For Kerr BHs of spin a, the following 
relation holds, 



4\/27r 



1 + VT^ 



1/2 



(28) 



and hence, a is determined from Ce and Min- We esti- 
mate the spin assuming that Eq. ([28|) holds even for the 
BH surrounded by the disk [44] . The BH spin determined 
by this method is referred to as a{2- 

In the third method, Cp/Ce is used. For Kerr BHs, it 
is calculated to give 



Cp 



-E{ay2f+), 



(29) 



where f-|_ = 1 
defined by 



Vl — and E{z) is an elliptic integral 



E[z) = / 
Jo 



z sin^ ede. 



(30) 



Thus, assuming that the same relation holds even for the 
BH surrounded by the disk, we may estimate a BH spin 
from Cp/Ce- We refer to this spin as a^. 

The values of afi-afa for TV = 36 are listed in Table 
III. We find that three values agree within a few % for 
Q = 3-5. The values of the spin depend weakly on the 
compactness of the NSs, and hence, we conclude that the 
spin parameter of the formed BH is ~ 0.56 ± 0.01, 0.48 ± 
0.01, and 0.42±0.01 for Q = 3, 4, and 5, respectively. For 
the smaller values of Q, the final BH spin is larger. The 
reason for this is that for the larger mass ratio, the total 
angular momentum of the system J at a given value of 
rriQ^l in the inspiral orbit is approximately proportional 
to rn^Q/iX + Q)^ (see also the initial condition in Table 
I). Thus, the binaries of smaller values of Q should form 
a BH of higher spin. 

For the case that Q < 2 and the disk mass is > O.OIM*, 
afi does not agree well with af2 and aa with the error 
size ~ 0.05. Because af2 and Ofs agree well, the error 
in flfi seems to be much larger than those of af2 and 
af3 . Then, the possible error sources are (i) underestima- 
tion of angular momentum that the disks possess and/or 
(ii) underestimation of angular momentum dissipated by 
gravitational waves. The possibility (ii) is not very likely, 
because for Q > 3 and for models M20.160 and M20.178, 
Ofi-flfa agree in a much better manner, indicating that 
gravitational waves are computed with a good accuracy. 
A possible reason that the angular momentum of the disk 
is underestimated is that many of the disk materials are 
located not in the finest grid domain in the AMR grid 



but in the second-fourth finest domains in which the grid 
resolution may not be high enough, and hence, the an- 
gular momentum is spuriously dissipated. This reason is 
also inferred from Fig. [7] (d) , which shows that the disk 
mass surrounding BH depends on the grid resolution. 



C. Gravitational waves 

1. Comparison with Taylor T4 waveform 

Figure[2](a)-(f) plot gravitational waveforms {+ mode) 
observed along the z axis as a function of the re- 
tarded time for models M20.145-M50.145, M20.160, and 
M30.178. {h denotes a gravitational wave amplitude.) 
The retarded time is approximately defined by 



t,,t = t- D-2Moln{D/Mo), 



(31) 



where D is a distance between the source and an observer. 
The gravitational waveforms shown here are obtained by 
performing the time integration of the Newman-Penrose 
quantity of I = \m\ = 2 mode. From the values of Dh/mo 
shown in Fig. [HI the amplitude of gravitational waves at 
a distance D is evaluated by 



.,.«2.4xl0-(^) 



100 Mpc\ / mo 



D 



bMr. 



(32) 



To validate the numerical waveforms presented here, 
we first compare the waveforms in the inspiral phase with 
those derived by the so-called Taylor-T4 formula for two 
point masses in quasicircular orbits. In the Taylor-T4 
formula, one calculates evolution of the angular veloc- 
ity, il, of the quasicircular orbits due to the gravita- 
tional radiation reaction up to the 3.5PN level beyond the 
quadrupole formula: The circular orbits at a given value 
of 17 are determined by the 3PN equations of motion ne- 
glecting gravitational radiation reaction, and then, one 
considers an adiabatic evolution of using the 3.5PN 
formula f or g ravitational radiation reaction (see, e.g., 
Refs. [i^ l46l | for a detailed description of various for- 
mulas based on the PN theory). Recent high-accuracy 
simulations for equal-mass (nonspinning or corotating) 
BH-BH binaries have proven that the Taylor-T4 formula 
provides their orbital evolution and gravitational wave- 
forms with a high accuracy at least up to about one orbit 
before the onset of the merger. Assuming that this holds 
for unequal-mass binaries, we calibrate our numerical re- 
sults by comparing them with the results by the Taylor- 
T4 formula. (Indeed, our numerical results indicate that 
the Taylor-T4 formula provides a good approximate so- 
lution for U,{t) even for the nonequal-mass binaries as 
shown below.) In the present work, the 3PN formula [i^l 
is employed for calculating the amplitude of gravitational 
waves in the Taylor-T4 formula. 
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TABLE III: Rest mass of material located outside apparent horizon (Mr>rAH)i mass estimated by energy-conservation 
law (MBH,f), BH mass estimated from equatorial circumferential radius (Ce/47r), irreducible mass of the BH (Afin; square 
root of area of apparent horizon in units of IGtt), ratio of polar circumferential radius (Cp) to equatorial one (Ce) of apparent 
horizon (i.e., Cp/Ce), and estimated spin parameters of the final state of the BH. an, an, and Sfa are computed from BH mass 
and angular momentum estimated by conservation laws, from Mirr and Ce of apparent horizon, and from Cp/Ce, respectively. 
All the values presented here are measured for the state obtained at the end of the simulations for = 36. Note that the 
parameters for the BH still vary with time at the end of the simulation for models M15.145, M20.145, and M20.145N because 
of gradual mass accretion. The error for Ce and Cp/Ce is ^ 0.1%, whereas that for an, at2, and afs is ^ 0.01 except for the 
case that the disk is formed for which the error of an would be ~ 0.05. The values for the final state of the BH depend very 
weakly on the grid resolution as far as A'^ > 30, but the mass of disk which presents only for Q < 2 systematically increases 
with N. The present results should be regarded as the lower-bound for the disk mass. 



Model 


Mr>rf,^/M» MBH,f/Mo 


Ce/47rMo 


Mirr /Mo 


Cp/Ce 


flfi 


af2 


af3 


M15.145 


0.023 


0.983 


0.981 


0.895 


0.872 


0.801 


0.747 


0.750 


M20.145 


0.010 


0.988 


0.984 


0.915 


0.898 


0.717 


0.684 


0.682 


M20.145N 


0.011 


0.988 


0.983 


0.916 


0.900 


0.721 


0.677 


0.676 


M20.160 


6 X 10"" 


0.988 


0.987 


0.919 


0.899 


0.694 


0.679 


0.680 


M20.178 


< lo-"" 


0.983 


0.983 


0.917 


0.904 


0.676 


0.672 


0.666 


M30.145 


< 10-5 


0.989 


0.985 


0.942 


0.934 


0.566 


0.559 


0.564 


M30.160 


< IQ-'' 


0.985 


0.983 


0.940 


0.935 


0.547 


0.560 


0.562 


M30.178 


< 10-5 


0.982 


0.981 


0.940 


0.937 


0.551 


0.550 


0.552 


M40.145 


< 10-5 


0.988 


0.986 


0.955 


0.953 


0.485 


0.482 


0.474 


M50.145 


< 10-5 


0.989 


0.986 


0.963 


0.964 


0.408 


0.419 


0.425 



TABLE IV: Several outputs of gravitational waves. Total radiated energy (AE) and angular momentum (A J) in units of initial 
ADM mass (Mo) and initial angular momentum (Jo), fraction of radiated energy for {I, \m\) = (2,2), (3,3), (4,4), and (2,1) 
modes, frequency of the fundamental quasinormal mode, kick velocity, and type of the gravitational waveform. The radiated 
energy for each mode is shown in unit of Mq by %. The origin of the error bar is primarily the numerical error associated 
with the finite grid resolution, and in part, the finite extraction radii of gravitational waves. Note that gravitational waves are 
extracted for several coordinate radii of 50-lOOAfo . 



Model 


AE/Mo (%) AJ/Jo (%) 


(2,2) 


(3,3) 


(4,4) 


(2,1) 


/qnmMbh Vkick (km 


) Type 


M15.145 


0.68 ±0.02 


14 ± 1 


0.66 ±0.02 


0.006 ± 0.002 


0.003 ± 0.001 


< 0.01 




< 5 


I 


M20.145 


0.87 ±0.02 


17.4 ±0.3 


0.85 ±0.01 


0.017 ±0.002 


0.005 ± 0.001 


0.003 ± 0.001 




21 ±5 


I 


M20.145N 


0.78 ±0.02 


15.0 ±0.2 


0.76 ±0.01 


0.014 ±0.001 


0.005 ± 0.001 


0.002 ± 0.001 




11 ±2 


I 


M20.160 


1.22 ±0.02 


22 ± 1 


1.19 ±0.02 


0.025 ± 0.005 


0.006 ± 0.002 


< 0.01 




62± 15 


II 


M20.178 


1.7 ±0.1 


25 ± 1 


1.6 ±0.1 


0.05 ±0.01 


0.01 ±0.005 


0.03 ±0.01 


0.087 


126 ± 20 


II 


M30.145 


1.3 ±0.1 


22.0 ±0.4 


1.2 ± 0.1 


0.08 ±0.01 


0.014 ±0.002 


0.03 ±0.02 


0.081 


98 ±9 


II 


M30.160 


1.6 ±0.1 


26 ± 1 


1.4 ±0.1 


0.09 ±0.02 


0.020 ± 0.003 


0.05 ± 0.02 


0.080 


153 ± 16 


HI 


M30.178 


1.8 ±0.1 


26 ± 1 


1.7±0.1 


0.12 ±0.01 


0.02 ±0.01 


0.04 ± 0.03 


0.080 


137 ± 43 


HI 


M40.145 


1.3 ±0.1 


23.5 ±0.5 


1.1 ± 0.1 


0.13 ±0.01 


0.025 ±0.005 


0.06 ±0.02 


0.077 


136 ± 12 


HI 


M50.145 


1.11 ±0.05 


24 ± 1 


0.85 ±0.02 


0.13 ±0.01 


0.04 ±0.01 


0.09 ±0.01 


0.074 


137 ±6 


HI 



More specifically, the comparison of the numerical 
waveforms and semianalytic ones derived by the Taylor- 
T4 formula is carried out in the following manner. First, 
we derive the orbital angular velocity as a function of 
time for a numerical result from ^4 by [isj 



n{t) = 



1 |*4(/ = ™ = 2)| 



dm 41 = m = 2) 



(33) 



where 'i/ 4(1 — m — 2) is the I = m = 2 mode of '^4. 
Then, we compare the numerical result for il(t) with the 
semianalytic one derived by the Taylor-T4 formula (see 



Fig, mi)) . When comparing two results, we have a degree 
of freedom for the time translation. Thus, first of all, by 
shifting the time axis of the Taylor-T4's result, we align 
the origin of the time. As shown in Fig. 1101 we can always 
shift the time axis appropriately to align two results for 

After the appropriate time translation, we compare the 
gravitational waveforms obtained by the numerical simu- 
lation and by the Taylor-T4 formula. Then, we still have 
a degree of freedom for choosing the wave phase. Thus, 
we iteratively change the phase of the Taylor-T4's wave- 
form until a good matching of two waveforms is achieved. 

In Fig. [51 the dot-dotted curve denote the resulting 
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FIG. 10: Orbital angular velocity computed from ^'4 as a function of the retarded time for models (a) M20.145, (b) M30.145, 
(c) M40.145, and (d) M50.145. tret denotes the retarded time and mo is the total mass. For all the panels, the results for two 
different grid resolutions with = 30 (dashed curve) and 36 (solid curve) and the results for M20.145b-M50.145b (dot-dashed 
curves) are shown together. The dot-dotted curve denotes the result derived by the Taylor- T4 formula. 



semianalytic waveform derived by the Taylor- T4 formula. 
It is found that the numerical waveforms agree with the 
Taylor-T4's results with a good accuracy, except for the 
early stage of the simulations (i.e., for the first a few wave 
cycles), during which the eccentricity of the binaries is 
not small. In particular, for the last several wave cycles 
(except for the orbit just before the merger), the numer- 
ical wave phases agree with those derived by the Taylor- 
T4 formula with in ~ 3% error. This indicates that the 
binaries computed in the present simulation are indeed in 
an approximately quasicircular orbit of small eccentric- 
ity at least for the last several inspiral orbits. Also, this 
indicates that the Taylor-T4 formula provides a good ap- 
proximate solution for f2(t) even for the nonequal-mass 
binaries, because the agreement systematically holds ir- 
respective of the mass ratio. 

Figure [TO] also shows good agreement between the nu- 
merical and Taylor-T4's results for n{t). The numerical 
results for two grid resolutions [TV — 36 (solid curve) 
and 30 (dashed curve)] are shown for illustrating that a 
convergence is approximately achieved. In addition, the 
results for models M20.145b-M50.145b are plotted for 
comparison. This figure shows that the evolution of il(t) 
for models M20.145-M50.145 agrees well with those pre- 
dicted by the Taylor-T4 formula for a long time duration 
irrespective of the mass ratio. Modulation in the evolu- 
tion of is Afi/fi ^ 6% for the last a few inspiral phase 
(irrespective of A'' = 30 or 36), and thus, the eccentricity, 
which is approximately estimated by 2A^l/3n, is ^ 4%. 
In the last several orbits, the eccentricity appears to be 



at most ^ 1%. 

Comparing the results of M20.145-M50.145 with those 
of M20.145b-M50.145b, we find that the modulation is 
by a factor of 2 larger with the initial condition com- 
puted in the [3^^ condition. This unfavored behavior is 
more significant for the larger values of Q. This demon- 
strates the advantage for using the initial condition com- 
puted in the 3PN-J condition. 



2. Classification of waveforms 

Figure [H] shows that gravitational waveforms in the 
merger and ringdown phases depend sensitively on the 
mass ratio and compactness of the NS. Comparison of the 
waveforms for models M20.145 and M20.160 [see Fig. [S] 
(a) and (e)], for which masses of the BH and NS are 
identical each other, illustrates a strong dependence of 
the merger waveforms on the NS compactness. The wave 
amplitude for model M20.145 decreases suddenly in the 
middle of the inspiral phase due to tidal disruption. By 
contrast, the wave amplitude for model M20.160 does 
not decrease as quickly as that for model M20.145 be- 
cause the tidal disruption docs not occur far outside the 
ISCO. 

Gravitational waveforms in the merger and ringdown 
phases for models M30.145 and M30.178 [see Fig. [9] (b) 
and (f), and Fig. fTTT d)] are also distinguishable because 
the amplitude in these phases is much larger for model 
M30.178. This is due to the fact that the NS for model 
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FIG. 11: Gravitational waves (the real part of ^'4) emitted in the merger and ringdown phases for models (a) M30.145, (b) 
M40.145, (c) M50.145, and (d) M30.145 (dashed curve), M30.160 (dotted curve), and M30.178 (solid curve). For (a)-(c), a 
fitting waveform given by Eq. (|34p is plotted together by the dot-dotted curve. 



M30.178 is not strongly affected by the tidal force of the 
companion BH even at the ISCO. By contrast, the tidal 
effects play an important role for deformation of the NS 
in close orbits for model M30.145. Because the NS is dis- 
rupted near the ISCO for this model, the wave amplitude 
in the merger and ringdown phases is significantly sup- 
pressed. These results illustrate that gravitational waves 
emitted in the merger and ringdown phases have a po- 
tential information about the compactness of the NS (see 
also Sec. HVDl) . 

As these comparisons clarify, there are three types 
of gravitational waveforms. For the case that the NS 
is tidally disrupted during the inspiral phase (e.g., for 
model M20.145), the wave amplitude quickly decreases 
at the tidal disruption, as the waveforms associated with 
the inspiral motion are suddenly shut off. Namely, the 
waveform is composed only of the inspiral waveform and 
subsequent sudden shut-off, and the merger and ring- 
down waveforms are essentially absent. We refer to this 
type as the type I. 

Even in this case, most of the material of the tidally 
disrupted NS falls into the companion BH (see Fig. [7]). 
During the tidal disruption and subsequent infalling into 
the BH, an orbital motion of the disrupted material and 
an oscillation of the BH may excite merger and ringdown 
gravitational waves (here "ringdown gravitational waves" 
imply gravitational waves associated with quasinormal 
modes of the BH). However, such waveforms are not seen. 
The likely reason for the absence of the merger waveform 
is that the NS is significantly elongated in a short time 
scale and its density quickly decreases, suppressing an ef- 



ficient excitation of gravitational waves. The reason for 
the absence of the ringdown waveform is that the mate- 
rial of the NS, which falls into the BH, does not have a 
compact configuration but have an elongated low-density 
configuration. In the case that such low-density diffuse 
matter incoherently falls into the BH, the excitation of 
the quasinormal modes is significantly suppressed due to 
the phase cancellation effect [48| . 

For models M30.145 and M20.160, mass shedding oc- 
curs before the binaries reach the ISCO. However, the 
sudden shut-off of the inspiral waveforms is not seen be- 
cause the tidal disruption and the subsequent spreading 
of the material do not occur during the inspiral phase. 
Rather, most part of the elongated NS falls into the BH 
before the tidal disruption is completed. In this case, 
the ringdown waveform is seen but the amplitude is low 
because the quasinormal mode is not excited efficiently. 
By contrast, just before the ringdown gravitational waves 
are emitted, gravitational waves are significantly excited 
by a matter motion around the BH. Thus, the merger 
gravitational waves are present. Namely, in these cases, 
the gravitational waveforms are composed primarily of 
the inspiral and merger ones. We refer to this type as 
the type II. 

For models M40.145, M50.145, and M30.178, tidal ef- 
fects to the NS do not play an important role. In this 
case, the gravitational waveform is composed of the in- 
spiral, merger, and ringdown waveforms, as in the merger 
of BH-BH binaries (e.g., [H, HI]). We refer to this type 
as the type HI. 
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3. Ringdown waveforms 

For models M30.145, M40.145, M50.145, M30.160, 
M20.178, and M30.178, ringdown gravitational waves as- 
sociated with quasinormal modes of the formed BH are 
excited in the final phase. These gravitational waves are 
approximately described by 



Ae-'/'" sin(27r/QNMt + 5), 



(34) 



where A and 5 are constants, and /qnm and t^i are the 
frequency and damping time scale of the fundamental 
quasinormal mode. A perturbation study predicts the 
frequency and damping time scale for a BH of mass A/bh 
and spin a as [49| 



/qnm-^bh 



0.16[1 

1-0.45 



0.63(1 -a)"-^] 



"■/qnm 



(35) 
(36) 



For models M30.145, M40.145, and M50.145, the BH spin 
is « 0.56, 0.48, and 0.42 as shown in Sec. HVBl (For 
models M30.160 and M30.178, the spin agrees approx- 
imately with that for model M30.145.) Thus, for each 
of these models, /qnmA/bh ~ 0.081, 0.077, and 0.074 
(/qnmMo « 0.082, 0.078, and 0.075), and nUfq^M « 
2.9, 2.7, and 2.6, respectively. 

Figure[TT](a)-(c) compares the numerical waveforms in 
the ringdown phase with the hypothetical analytic wave- 
forms for models M30.145-M50.145. This shows that 
the numerical waveforms are fitted by the analytic one 
([M)) fairly well, and that /qnm and td computed from 
the data of the BH geometry agree approximately with 
those computed from gravitational waveforms. However, 
the numerical waveforms do not agree completely with 
the hypothetical ones. This disagreement is reasonable 
because for these models, the material of the NS does 
not simultaneously fall into the BH at the merger. Thus, 
gravitational waves are emitted both by a motion of the 
material moving in the vicinity of the BH and by the 
quasinormal-mode oscillation of the BH. In addition, the 
quasinormal modes are not simultaneously excited be- 
cause the material does not simultaneously fall into the 
BH, and the resulting waveforms may be composed of 
many ringdown waveforms, as well as of the waveforms 
excited by a material moving around the BH. Further- 
more, the system is not completely in a vacuum nor in 
a stationary state, and hence, the numerical waveforms 
may not be fitted precisely by the analytic results derived 
in the ideal assumption. 

The wavelength of gravitational waves emitted in the 
merger phase (around the time when the peak amplitude 
is reached) is slightly longer than that of the quasinormal 
mode (i.e., the frequency is lower than /qnm)- This in- 
dicates that these gravitational waves are not emitted by 
an oscillation of the BH, but they are likely to be emitted 
primarily by a motion of the material which moves in the 
vicinity of the BH. In the merger phase, the amplitude 
gradually decreases after the peak is reached. The reason 



for this behavior is that the material is elongated by the 
tidal effect of the BH during the infalling; i.e., the rea- 
son is not the damping associated with the quasinormal 
mode oscillation. The amplitude emitted in the merger 
phase is much larger than that emitted in the ringdown 
phase, although the characteristic frequencies of these 
two types of gravitational waves are not very different. 
In Sec. IIVDI we find that the Fourier spectrum has a 
plateau in a high frequency region /mp ~ 0.06-0.08 for 
the case that Q = 3-5. This plateau is primarily gener- 
ated by gravitational waves emitted in the merger phase 
not by the quasinormal mode oscillation. 

Figure [TT] (d) compares the waveforms emitted in 
the merger and ringdown phases for models M30.145, 
M30.160, and M30.178. For these models, the masses 
of the BH and NS are identical each other, but the com- 
pactness of the NSs is different. This figure shows that 
the amplitude is larger for the model with the larger value 
of C. This is reasonable because more compact NSs are 
less subject to the tidal effects by the companion BH, and 
hence, the material of the NS falls into the BH in more 
simultaneous manner, resulting in a coherent excitation 
of gravitational waves. 

The difference in the amplitude of gravitational waves 
emitted in the merger and ringdown phases is reflected 
also in the noticeable difference of energy and angular 
momentum carried by gravitational waves. As shown 
in Table IV, for example, the total energy radiated for 
models M30.160 and M30.178 is by ~ 35% and - 50% 
larger than that for model M30.145, respectively. 



D. Gravitational wave spectrum 

To determine the effective amplitude of gravitational 
waves for a given frequency, the Fourier spectrum of grav- 
itational waves of I ~ |m| = 2 modes are computed. In 
this paper, as the Fourier spectrum, we define 



Kf) 



\h+{fW + \h^{fW 



(37) 



where / is the frequency, h+{f) and (/) are the Fourier 
transformation of the plus and cross modes of gravita- 
tional waves observed along the z axis; 



(38) 
(39) 



Then, the most optimistic effective amplitude of gravita- 
tional waves for a given frequency is defined by fh(f). 

In the numerical simulations, BH-NS binaries of a fi- 
nite orbital separation are prepared as the initial con- 
dition and the inspiral phase is computed for a finite 
duration. Consequently, the Fourier spectrum for the 
low- frequency side (for / ^ f^o /'"') becomes absent if we 
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FIG. 12: (a)-(d) The spectrum of gravitational waves fh{f)D/mo (solid curves) for models (a) M20.145 and M20.160, (b) 
M30.145, (c) M40.145, and (d) M50.145, respectively. The dashed and long-dashed curves denote the relation according to Eq. 
(|40[1 and spectra of gravitational waves computed in the Taylor- T4 formula. The dotted curves denote the results of fitting. The 
upper horizontal and right vertical axes show the value in a hypothetical value of Mns = l.SSM© and D = 100 Mpc. The arrow 
indicates the frequency of the fundamental quasinormal mode calculated by Eq. (|36|) . (e) The same as (a) but for gathering the 
spectra for models M20.145-M50.145 (long-dashed, dotted, dashed, and solid curves). To clarify the qualitative feature of the 
spectra, a smoothing procedure is applied for the numerical data. The three lines in the upper left side denote the predicted 
frequency at which mass shedding of the NS occurs due to tidal field of the companion BH for models M20.145-M40.145 {Q 
denotes the mass ratio), (f) The same as (b) but for models M30.145 (dot-dotted curve), M30.160 (dotted curve), and M30.178 
(solid curve). A smoothing procedure is also applied for the numerical data. 
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naively perform the Fourier transformation for the nu- 
merical data. To compensate the Fourier spectrum for 
the low-frequency side, we combine a hypothetical wave- 
form computed by the Taylor-T4 formula, as often done 
(e.g., Refs. [3, S^)- To do so, wc match the numeri- 
cal waveforms with those by the Taylor-T4 formula at a 
time when the corresponding binary orbit in the numer- 
ical simulation is approximately in a quasicircular state. 
As shown in Fig. [9l two waveforms match well for a wide 
range of time, so the resulting Fourier spectrum depends 
only weakly on the chosen matching time. 

Figure [T^] plots fh{f)D/mo for various models. In 
Fig. [12] (a)-(d), we also plot the Fourier spectra of a 
gravitational waveform derived in the Taylor-T4 formula 
(long-dashed curve) and by the Newtonian waveform 
(dashed curve) as (e.g., (s^l) 



fKf) 



5 mo 



)l/2 
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(Trmo/) 



-1/6 



(40) 



Here, "Newtonian" implies that the orbital motion is cal- 
culated in the Newtonian plus 2.5 PN equations of mo- 
tion, and the gravitational waveform is computed by the 
quadrupole formula. 

As Figure [T^] indicates, the spectra of gravitational 
waves emitted in the inspiral, merger, and ringdown 
phases are smoothly connected. Nevertheless, we still see 
modulations of small amplitude for 0.01 < /toq ^ 0.02, 
which are likely to be caused by slight modulations in 
the amplitude and/or phase of numerical gravitational 
waveforms. 

In the upper-horizontal and right-vertical axis, we plot 
the frequency and averaged effective amplitude for hypo- 
thetical values Mns = 1.35M0 and D = 100 Mpc. Here, 
the averaged effective amplitude is defined by the aver- 
age of fh(f) over the source direction and the direction 
of the binary orbital plane as 



/leff = OAfhif) 



= 9.6 X 10 
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Figure fT2l shows that the spectrum shape has the fol- 
lowing universal qualitative feature: Irrespective of the 
values of C and Q, fh{f) (hereafter referred to as the 
spectrum amplitude) decreases with the increase of / 
and above a "cut-off" frequency, /cut, it decreases ex- 
ponentially. For / — > 0, fh{f) is universally propor- 
tional to /^^/®, and / < /cut, fh{f) is always written as 
oc f~^^^F{f), where F{f) is a slowly varying function of 
/ and the value is ^ 1. However, the detailed spectrum 
shape depends on the values of C and Q, and Fig. [T2] ex- 
hibits a wide variety of the possible spectrum shape as 
explained in the following. 

For the case that the NS is tidally disrupted outside 
the ISCO, type I gravitational waves are emitted (e.g., for 



model M20.145). In this case, the spectrum amplitude 
monotonically decreases, and above a cut-off frequency, 
it exponentially decreases with the increase of /. The 
cut-off frequency for model M20.145 is /cut ~ 0.04/too. 
For / < 0.03/too ~ 0.8/cut, Fif) may be approximately 
fitted by the 3PN formula for the amplitude [13], and 
thus, the spectrum may be described, e.g., by the follow- 
ing way 



Hf) = /^3PN(/)e 



-(///c 



(42) 



where /13PN is the amplitude of gravitational waves de- 
rived in the 3PN theory (see Eq. (79) of Ref. [i^), and cr 
is a constant of order unity. The dotted curve in Fig. [12] 
(a) denotes the result of the fitting for /cut = 0.038/too 
and a = 2.2, and we see that the fitting works well. 
For model M15.145, the fitting also works well and gives 
/cut = 0.030/too and a — 2.2. Because the NS is tidally 
disrupted farther outside the ISCO, the value of /cut is 
smaller for model M15.145. 

The latest high-precision numerical study for the 
quasiequilibrium states of the BH-NS binaries shows that 
the mass shedding of the NSs occurs if the following con- 
dition is satisfied 7] 
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where the value of C is ~ 0.270 for the F = 2 polytropic 
EOS. (Note that this is the relation which holds only for 
the NS of the irrotational velocity field and for the BH 
of no spin.) From this relation, we expect that the mass 
shedding sets in at a frequency /ms = ^^ms/tt- We plot 
this expected frequency (0.0174, 0.0219, and 0.0265/too) 
for models M20.145-M40.145 together with the Fourier 
spectra for models M20.145-M50.145 in Fig.[n](e). 

Figure [12] (a) and (e) show that any special feature is 
not seen at / = /ms in the spectrum even for the case 
that the NS is tidally disrupted before the binary reaches 
the ISCO, e.g., for model M20.145 (a small bump seen for 
/ < /ms is due to modulation of the wave amplitude and 
wave frequency in the numerical data). Also, it is found 
that /cut is much larger than /ms, a-s already pointed out 
in Ref. [l4|. This implies that at / = /ms, the mass- 
shedding sets in for the NS, but the NS is not tidally 
disrupted at such a low frequency. Even for the closer 
orbits of / > /ms, the NS behaves as a self-gravitating 
star although it transfers a small amount of mass to the 
companion BH, and gravitational waves from this BH-NS 
binary are characterized basically by the inspiral wave- 
forms. However, because the tidal effect becomes more 
important for the smaller orbital separation, the tidal 
deformation of the NS is enhanced more with time, and 
eventually, the tidal disruption occurs. At the tidal dis- 
ruption, the gravitational wave amplitude should quickly 
decrease, so it is natural to identify the frequency of the 
last inspiral wave as /cut in this case. 

Here, we should emphasize that the compactness (or 
radius ) of the NSs is reflected in /cut not in /ms • In the 
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existing idea, one naively assumes that the tidal disrup- 
tion occurs at / = /ms and discusses that the radius 
of the NS (i.e., the EOS of the NS) may be constrained 
by identifying the values of /ms (e.g. Ref. [5?|). How- 
ever, our present result indicates that this possibility is 
unlikely. To constrain the EOS of the NS from gravita- 
tional waves detected, we have to theoretically determine 
/cut for a variety of the EOSs and masses of the BH and 
NS. 

For models M20.160 and M30.145 [see Fig. [11(a) and 
(b)], the spectrum shape is similar to that for model 
M20.145, but the cut-off frequency is higher as /cut ~ 
0.06/too (see below). This cut-off frequency is still lower 
than the frequency of the fundamental quasinormal mode 
(/qnm; see the arrow of Fig. [H (a) and (b)). This in- 
dicates that the tidal disruption occurs at an orbit cor- 
responding to f ~ /cut- This cut-off frequency, which is 
higher than that for model M20.145, reflects the fact that 
the tidal disruption occurs at a closer orbit than that for 
model M20.145 due to the larger compactness or mass 
ratio of the system. The cut-off frequency appears to be 
much higher than the frequency of gravitational waves 
emitted at the ISCO, which is f\so ~ 0.03/mo. This im- 
plies that gravitational waves of fiso ^ / ^ /cut are not 
emitted by the inspiral motion of the binary but by a 
motion in the merger phase. (This is the reason why we 
do not classify the waveforms for models M20.160 and 
M30.145 into type I but into type IL) In this case, it is 
not appropriate to fit the spectrum by Eq. (|42p but by a 
different function (see below). 

Another interesting feature found in the spectrum of 
model M30.145 is that no peak associated with the quasi- 
normal mode appears at / = /qnm (see arrow in Fig. [H 
(b)). The reason is simply that the amplitude of the ring- 
down waveform is much smaller than that of the merger 
waveform, as pointed out in Sec. IIV CI 

For models M40.145, M50.145, and M30.178 (see 
Fig. [H (c), (d), and (f)), the gravitational waveforms 
are type III. In these cases, the spectrum amplitude also 
steeply decreases above a cut-off frequency, but the fea- 
ture of the spectrum shape is qualitatively different from 
those of models M20.145, M30.145, and M20.160, be- 
cause the gradual decrease of the spectrum amplitude 
continues approximately up to the frequency of the quasi- 
normal mode (see the arrows in these figures). Namely, 
/cut is approximately equal to /qnm, and thus, the cut- 
off frequency does not indicate that the tidal disruption 
occurs at an orbit of / = /cut- In other words, the signal 
of the tidal disruption is absent in these cases. 

Another difference is seen in the spectrum shape 
for a frequency slightly smaller than /cut- For mod- 
els M20.145, M30.145, and M20.160, a spectral index, 
n = -d\n{fh{f))/d\n{f), for / < /cut monotonically 
increases with /. For models M40.145, M50.145, and 
M30.178, by contrast, the value of n slightly decreases 
with /, and a "plateau" appears. This spectrum shape 
is similar to that for the merger of unequal-mass BH- 
BH binaries [5l| . This is reasonable because in this case. 



tidal effects do not play an important role during the in- 
spiral and merger phases, and the merger process should 
be qualitatively the same as that in BH-BH binaries. 

In Fig.[TKf), we compare the spectrum shape for mod- 
els M30.145, M30.160, and M30.178, for which the masses 
of the BH and NS are identical whereas the compactness 
of the NS is different. This figure illustrates that for 
the larger compactness of the NS, the cut-off frequency 
is higher. Also, the width of the plateau is larger for 
the more compact models (M30.160 and M30.178). As 
discussed above, these differences reflect the difference 
in the evolution process of the late inspiral and merger 
phases. Because the spectrum shape near / = /cut is sig- 
nificantly different among three models, observing such 
high-frequency gravitational waves will play a special role 
in constraining the compactness (i.e., radius) of the NS. 

To systematically identify the cut-off frequency, the fit- 
ting of the spectrum with a specific function is useful, as 
illustrated for model M20.I45. The spectra associated 
with the inspiral phase for models M30.145-M50.145, 
M20.160, M20.178, M30.160, and M30.178 are also fitted 
by the 3PN amplitude, hsp^if), for / < /iso ~ 0.03/mo. 
However, this should not be the case for /iso ^ / ^ /cut 
because gravitational waves for this high-frequency com- 
ponent are not emitted by the inspiral motion, but by a 
motion of the material associated with the merger and 
infalling into the BH. Hence, the spectrum should not be 
fitted by Eq. dH]) but by 

Hf) = /i3PN(/)e-(^//'"='" + /l„.crgcr(/), (44) 

where /imcrgcr(/) denotes the spectrum of gravitational 
waves emitted in the merger phase, and /ins is a fre- 
quency of 0.01-0. 03/mo. In this paper, we fix cr = 3.5 to 
reduce the total number of free parameters, and choose 
the following function for /imorgci (/), 

Wr(/) - [1 - e-(//^--)^], (45) 

where A and (Tcut are free parameters. We add a free pa- 
rameter /ins2 for which the value is close to /ins, because 
the fitting is achieved in a better manner with it. 

The dotted curves in Fig. [12] (b)-(d) denote the re- 
sults for (/ins"io, /ins2mo. A, fcutmQ, (Tcut) =(0.014, 0.014, 
0.130, 0.063, 2.90) for model M30.145, (0.016, 0.016, 
0.103, 0.079, 4.60) for model M40.145, and (0.019, 0.020, 
0.090, 0.087, 3.70) for model M50.145. As Fig. [H shows, 
the spectrum is well fitted by this simple fitting function. 

The fitting for model M50.145 (and M30.178) is not 
as excellent as those for other models. In these cases, 
the inclination of the plateau is rather flat in the high- 
frequency region. This feature is universal for the gravi- 
tational waveforms emitted in the merger of BH-BH bi- 
naries. Thus, in such a case, another fitting function 
proposed, e.g. in Ref. [5l|, may be better. The fitting 
procedure here is, however, still robust in extracting im- 
portant physical information, as discussed in the follow- 
ing. 
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FIG. 13: /cut as a function of C for various values of Q. The 
dashed and dotted lines denote the value of /qnm for Q — 
2 and 3, respectively. For Q = 4 and 5, /qnm is slightly 
smaller than that for Q — 3, and ~ 0.078/mo and 0.075/mo, 
respectively. 



the optimistic direction of the source with respect to the 
plane of the detector's arm and of its binary orbital plane, 
the effective amplitude is at most « 2.5x 10^^^. The de- 
signed sensitivity of the advanced LIGO is < 3 x 10"^^ at 
/ > 1 kHz [53] . This implies that it will not be able to de- 
tect gravitational waves during the merger and ringdown 
phases by the laserinterferometric detectors of standard 
design for D > 100 Mpc. To detect gravitational waves at 
such high frequency, the detectors composed of special in- 
strument (e.g., resonant-side band extraction), which has 
a high sensitivity for high-frequency gravitational waves, 
is necessary. As mentioned above, gravitational waves for 
/ ^ /cut carry important information of the NS radius. 
It is strongly desired that a detector which is sensitive for 
the frequency of 1 kHz ^ / ^ 3 kHz will be developed in 
the future. 



F. Kick velocity 



E. Cut-off frequency and determining the 
compactness of the NS 

In the fitting procedure described in the previous sub- 
sections, the most important output is the cut-off fre- 
quency, /cut, in which information about the compact- 
ness (or radius) of the NS may be reflected. For the 
case that the NS is swallowed by the BH with no tidal 
disruption (e.g., Q ^ 4), it approximately indicates the 
frequency of the quasinormal mode of the formed BH, 
and thus, it will not be possible to extract the compact- 
ness of the NSs from the gravitational-wave signal. By 
contrast, for the case that the NS is tidally disrupted be- 
fore it is swallowed by the BH (i.e., for the small values 
of Q and C), /cut is much smaller than /qnm and in- 
dicates a characteristic frequency of gravitational waves 
emitted at the tidal disruption (see Fig. [T5)) . We find 
that for models M20.145, M20.160, M20.178, M30.145, 
and M30.160, /cutmo « 0.038, 0.056, 0.070, 0.063, and 
0.083. For Q — 2, the value of /cut is smaller than the 
frequency of the quasinormal mode, /qnm ~ 0.09/mo for 
a wide value of C. This indicates that we will be able to 
constrain the compactness of the NS if we detect gravita- 
tional waves in the merger phase for a binary composed 
of a low-mass BH and an NS with Q ^ 2. For model 
M30.160, /cut is approximately equal to /qnm, whereas 
/cut < /qnm for mode M30.145. Thus, if the compact- 
ness of the NS is fairly small C < 0.16, we will have a 
chance that the detection of gravitational waves from a 
binary of Q w 3 constrains the compactness of the NS. 
By contrast for Q > 4, gravitational waves are not likely 
to have robust information about the compactness of the 
NSs. 

Figure [12] (a)-(d) and (f) show that the cut-off fre- 
quency is / = /cut ^ 1.4-2.2 kHz and the effective ampli- 
tude at / = /cut is universally 1 x 10~^^ for hypothet- 
ical values D = 100 Mpc and Mns = 1.351/©. Even for 
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FIG. 14: Kick velocity as a function of u = Q/{1 + Q)^ for 
different compactness of the NSs. The dot-dotted curve de- 
notes the fitting formula derived in Ref. [iOl for the merger of 
nonspinning BH-BH binaries. 

Kick velocity induced by anisotropic gravitational ra- 
diation is estimated from the total linear momentum car- 
ried by gravitational waves. Figure [TH plots the kick ve- 
locity as a function of the ratio of reduced mass to total 
mass for different compactness of the NSs (see also Table 
IV). 

We find that the kick velocity depends strongly on the 
mass ratio and compactness of the NS. More specifically, 
it depends strongly on whether the NS is tidally disrupted 
or not. If the tidal disruption does not occur, the NS 
simply falls into the companion BH in a basic picture 
(although tidal deformation and strong elongation of the 
NS still occur just before it falls into the BH). In such a 
case, the merger process is similar to that in the merger 
of nonspinning BH-BH binaries, and the kick velocity 
should agree approximately with the results for those 
systems. To confirm this fact, we compare our numer- 
ical results with the fitting formula derived in Ref. [i^l 
(see dot-dotted curve in Fig. [T?)l . Figure [U shows that 
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our results for models M40.145, M50.145, M30.160, and 
M30.178 are in a reasonable agreement with the fitting 
formula of Ref. [13]; the kick velocity is - 100-200 km/s. 
This indicates that our numerical results are reliable. 

In the case that the NS is tidally disrupted, the kick 
velocity is significantly suppressed. The primary reason 
is that the kick is excited most efficiently when two stars 
merge; more specifically, when the amplitude of gravita- 
tional waves becomes maximum. The maximum ampli- 
tude of inspiral gravitational waves is higher for the case 
that two stars have more compact orbits (i.e., the NS 
is more compact). Also, the amplitude of gravitational 
waves in the merger phase becomes higher when an NS 
falls into a companion BH without tidal deformation. For 
the case that the tidal disruption happens, close inspiral 
orbits are prohibited, and as a result, the enhancement 
of the gravitational-wave amplitude is suppressed as we 
showed in Sec. IIV CI Due to this fact, the kick velocity 
is smaller for the smaller values of compactness and mass 
ratio. 



V. SUMMARY 

New general relativistic simulation for the inspiral, 
merger, and ringdown of BH-NS binaries are performed 
by a new AMR code SACRA. This paper presents the nu- 
merical results of a longterm simulation for a variety of 
mass ratio and compactness of the NS for the first time. 
For the simulation, we adopt initial conditions for which 
the unphysical eccentricity is much smaller than that in 
our previous study [l3|, and also the initial orbital 
separation is much larger. In the present results, the 
binaries spend in the inspiral phase for 4-7 orbits, and 
in the last several orbits, the eccentricity is sufficiently 
suppressed to be ~ 0.01 and the orbit becomes approxi- 
mately quasicircular. This is reflected in successful com- 
putation of gravitational waveforms for the inspiral phase 
which agree well with those predicted by the FN theory. 
Because the moving-puncture approach is adopted, the 
merger and subsequent evolution of the BH with the sur- 
rounding material are simulated for a long time until the 
system relaxes to a quasisteady state. As a result, accu- 
rate gravitational waveforms for the late inspiral, merger, 
and ringdown phases are derived successfully. 

By the present self-consistent simulation, we reconfirm 
the finding in our previous paper that the merger 
process depends sensitively on the mass ratio, Q, and 
compactness of the NS, C. Only for the case that both Q 
and C are sufficiently small, the NS is tidally disrupted 
by the companion BH. For example, for the chosen EOS 
in this paper (polytropic EOS with F = 2), the NS of 
C — 0.145 is tidally disrupted only for Q < 3 before 
the binary reaches the ISCO. By analyzing the spectrum 
of gravitational waves, we confirm that even in the case 
that the NS satisfies the condition for the onset of the 
mass shedding [see Eq. (|43|)]. the tidal disruption does 
not occur immediately. Rather, the NS behaves as a 



self-gravitating star for a while during the mass-shedding 
phase. Therefore, the tidal disruption occurs for the more 
restricted case than the mass shedding does. This fact 
also implies that for determining the condition of tidal 
disruption, numerical-relativity simulation is the unique 
approach. 

The parameter space for the formation of a long-lived 
disk surrounding the BH is even more restricted. For 
C = 0.145, the disk of lifetime ^ 10 ms is formed only 
for Q 2, and for C = 0.160, the disk mass is at most 
10~^M* even for Q = 2. The mass of the formed disk 
is > 0.01-0.02A/* for Q = 2-1.5 even for relatively less 
compact NSs with C — 0.145, and thus, the disk is not 
very massive, even if it is formed. This conclusion is 
consistent with the results of Ref. [l7| . 

At present, the precise EOS of high-density nuclear 
matter is unknown [43 |. However, many EOSs predict 
that the radius of NSs of canonical mass 1.25-1.45Mo 
is smaller than ~ 12 km (e.g., Ref. [13]). Namely, the 
compactness of the NS is larger than 0.155. If the EOS 
is really stiff and the NS is compact as the nuclear theory 
predicts, the disk will not be formed after the merger be- 
tween nonspinning BH and NS, as pointed out by Miller 
[ist . (However, this will not be the case if the BH has 
the spin of substantial magnitude.) 

We find that gravitational waveforms from BH-NS 
binaries are roughly classified into the following three 
types, (i) When the NS is tidally disrupted during the in- 
spiral phase, the gravitational waveform is characterized 
by the inspiral waveform and subsequent sudden shut- 
down. In this case, the amplitude of the Fourier spec- 
trum monotonically decreases with the increase of the 
frequency, and at a cut-off frequency /cut, which is the 
frequency of gravitational waves when the NS is tidally 
disrupted (not equal to /ms)i the spectrum amplitude de- 
creases exponentially [see Eqs. (|42p and (|44)l ]. We refer 
to this waveform as type I in this paper, (ii) In the case 
that the NS is tidally disrupted but most of the material 
falls into the companion BH before the tidal disruption 
is completed, the gravitational waveform is characterized 
by the inspiral and merger waveforms. In this case, ring- 
down gravitational waves are excited in the final phase, 
but its amplitude is much smaller than those of late in- 
spiral and merger gravitational waves. As a result, the 
shape of the Fourier spectrum is similar to that of type 
I but the cut-off frequency does not correspond to the 
frequency of the last inspiral orbit nor the frequency of 
the quasinormal mode. We refer to this waveform as 
type II. (iii) When the NS is not tidally disrupted before 
it is swallowed by the BH, the gravitational waveform 
is characterized by the inspiral, merger, and ringdown 
waveforms. In this case, the amplitude of the Fourier 
spectrum monotonically decreases with the increase of 
the frequency in the inspiral phase as in the cases (i) and 
(ii). However, in the late inspiral and merger phases, the 
gravitational- wave amplitude increases and, as a result, 
a plateau appears in the Fourier spectrum of / ^ /cut- 
Then the spectrum amplitude exponentially decreases for 
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/ > /cut [see Eq. (04])]. In this case, /cut is approximately 
equal to the frequency of the fundamental quasinormal 
mode of the formed BH, /qnm, and does not have any in- 
formation about tidal disruption (and thus compactness 
of the NS). We refer to this waveform as type III. 

We fit the spectrum of gravitational waves with hypo- 
thetical functions of a small number of parameters. We 
find that the fitting works well irrespective of the values 
of Q and C. By this fitting, we systematically determine 
the value of the cut-off frequency, /cut and confirm the 
following fact: For the case that the NS is not tidally 
disrupted, the value of /cut is approximately equal to the 
value of /qnm, as mentioned above. By contrast, for the 
case that the tidal disruption occurs, the value of /cut is 
smaller than /qnm, and it depends strongly on Q and 
C. The tidal disruption occurs for a wide range of C if Q 
is smaller than ^ 2 . Thus, if gravitational waves from 
binaries composed of a low-mass BH and NS at tidal dis- 
ruption are observed, it may be possible to constrain the 
compactness of the NS, and as a result, the EOS of the 
NS. 

We also reconfirm that the frequency at the onset of 
the mass shedding of the NSs is not reflected in the spec- 
trum in an outstanding manner. Namely, the compact- 
ness (or radius) of the NSs is reflected in /cut not in /ms- 
In the existing idea, one naively assumes that tidal dis- 
ruption occurs at / = /ms and discusses that the radius 
of the NS (i.e., the EOS of the NS) may be constrained 
by identifying the values of /ms (e.g. Ref. fs^l). How- 
ever, our present result indicates that this possibility is 
unlikely. To constrain the EOS of the NS from gravita- 
tional waves detected, we have to determine /cut for a 
variety of the EOSs and masses of the BH and NS, that 
can be done only by numerical-relativity simulation. 

Kick velocity induced by asymmetric gravitational 
wave emission is also computed. When tidal disruption 
does not occur (e.g., for (5 = 5), our result agrees approx- 
imately with that derived for the merger of nonspinning 
BH-BH binaries; the kick velocity is ~ 100-200 km/s. 
For the case that the tidal disruption occurs, the kick ve- 
locity is significantly suppressed, and it is much smaller 
than 100 km/s. This is due to the fact that gravitational 
wave amplitude in the last inspiral and merger phases is 



suppressed. 

As shown in this paper, gravitational waveforms in the 
merger and ringdown phases depend strongly on the com- 
pactness of the NS and mass ratio of the binary. This re- 
sults primarily from the fact that the degree of the tidal 
deformation and condition for the onset of the mass shed- 
ding and tidal disruption depend strongly on these two 
parameters. Gravitational waveforms should also depend 
on the EOS of the NS and on the spin of the BH, as illus- 
trated, e.g. in Ref. [5^1 ■ Thus, in the subsequent work, 
we will systematically perform the simulation for a vari- 
ety of the EOSs of the NS and BH spin, to clarify the 
dependence of the gravitational waveform and its spec- 
trum on these additional parameters. 

Note added in proof. After this paper was completed, 
we noticed that Illinois group posted a preprint which 
presented their latest results on the inspiral and merger 
of BH-NS binaries [5^. Gravitational waveforms they 
presented are qualitatively in good agreement with ours. 
However, in their results, the mass of disk surrounding 
a BH, which is formed after tidal disruption of the NS, 
is much larger than our results. For example, for Q — 3 
and C — 0.145 (with no spin for the BH), the disk is 
not formed in our results irrespective of the grid resolu- 
tion, but in their result, the disk mass is of order O.OlAf*. 
Currently, we do not understand the reason for this dis- 
crepancy. 
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